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ABSTRACT 



The dissertation constitutes the first detailed theoretical and experimental investigation 
of a thermoacoustic prime mover with periodic boundary conditions. There are five 
significant aspects to this research: (1) using DeltaE to analyze an annular prime mover, 

(2) developing an entirely new analysis program using MATLAB, (3) designing, building, 
and experimentally investigating a single stack, annular prime mover, (4) experimentally 
investigating a constricted, single stack annular prime mover, (5) predicting the 
performance of a two stack annular prime mover. The major conclusions are: (1) A single 
stack annular prime mover will not reach onset because the eigenmodes of the system do 
not support thermoacoustic growth. (2) A constricted annular prime mover will reach onset 
because the constriction forces dominating boundary conditions that alter the eigenmodes. 

(3) A two stack prime mover is predicted to reach onset because one of the eigenmodes of 
the symmetric system does support thermoacoustics. 
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I. INTRODUCTION 



The subject of this dissertation is the investigation of thermoacoustic prime movers in 
an annular geometry. What sets this research apart from previous work in prime movers is 
the nature of the boundary conditions. To our knowledge this work is the first thorough 
investigation of thermoacoustic prime movers with periodic boundary conditions. The 
primary conclusion drawn from this work is that a full understanding of the eigenmodes of 
the system are required to design a functioning annular prime mover. 

Thermoacoustic heat transport is a process through which an acoustic field generates, 
or is generated from, a flow of heat. Thermoacoustic engines are of two types: heat pumps 
and prime movers. A thermoacoustic heat pump utilizes a standing wave to transport heat 
along the boundary of a plate situated in the standing wave. The acoustically generated heat 
flow produces a temperature gradient across the stack. In other words, acoustic energy is 
converted into stored thermal energy. In contrast, a thermoacoustic prime mover produces 
acoustic work, by accepting heat from a high-temperature source and transferring it to a 
low-temperature sink. In this dissertation a prime mover in an annular resonator is 
investigated both theoretically and experimentally. 

To illustrate how a prime mover operates, two conventional prime mover 
configurations are described. One type (as is shown in Fig. 1.1) is a rigid-rigid standing 
wave tube which contains a stack of plates called the prime mover stack (or, simply, the 
stack), which is in thermal contact with two heat exchangers. In operation, a temperature 
difference is imposed between the two heat exchangers, resulting in a temperature gradient 
along the stack. When the temperature gradient in the stack is sufficiently large, the gas in 
the resonator oscillates spontaneously at a certain frequency, with pressure antinodes at the 
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closed ends. The inception of spontaneous oscillation is known as onset. Once onset is 
reached, the acoustic amplitude in the resonator grows rapidly, typically reaching several 
percent of the mean gas pressure. 




Figure 1.1 A typical rigid-rigid prime mover configuration. 

The second prime mover configuration, called the Hofler tube, is a rigid-open 
resonator (as is shown in Fig. 1.2). In 1983, Hofler designed and built this simple 
thermoacoustic oscillator to demonstrate some thermoacoustic phenomena to his doctoral 
committee at the University of California, San Diego. The sound radiated from the open 
end of the pipe is impressively loud, about 100 dB re 20 fiPa a meter away (Wheatley et 
al., 1985; Swift, 1988). 
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Figure 1 .2 The Hofler tube, an example of a rigid-open prime mover. 

The condition necessary for a prime mover to reach onset is that the amount of stored 
thermal energy converted to acoustic energy must exceed the total amount of acoustic 
energy dissipated by thermal-viscous losses in the prime mover. This ability to reach onset 
is determined by the geometry of the prime mover, in particular, the position of the stack 
relative to the nearest pressure antinode. The hot end of the stack is closer to a pressure 
antinode than is the ambient end. In the two examples of conventional prime movers 
described above, the node and antinode positions are fixed by incorporating well-defined 
acoustic boundary conditions into the system. This is an important common feature for all 
conventional prime movers. The “built-in” boundary conditions also serve as convenient 
starting points for computational analysis of the performance of a prime mover. Because 
the relative positions of the dominating boundaries and the stack are fixed, the stack 
position is optimized for only one acoustic mode. Also, the optimal stack position for 
onset is not necessarily the optimal position for high amplitude performance. One of the 
initial motivations for studying annular prime movers was to see if it would self-optimize 
the stack location relative to the acoustic field. 
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Unlike a conventional prime mover, an annular prime mover (as is shown in Fig. 1.3) 
does not have the typical dominating boundary condition to force a pressure or velocity 
node at any particular position. Design and analysis of a functional annular prime mover 
are thus made more difficult. 




Figure 1.3 An annular prime mover configuration. 

In a uniform cross-section annular resonator, the fundamental longitudinal mode 
corresponds to the circumference being approximately equal to one acoustic wavelength. 
There are two degenerate orthogonal modes that satisfy this condition. However, these 
modes have no preferred orientation. The presence of a nonuniformity (for example, the 
stack) breaks the degeneracy resulting in a frequency-splitting. The low frequency mode 
will have a pressure node at the stack, and the higher frequency mode will have a velocity 
node at the stack. (In what follows, the word “frequency” will be dropped and the modes 
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referred to, simply as “high” and “low”.) Neither of these modes supports 
thermoacoustics. 

Although in retrospect the answer now seems obvious, prior to a thorough study of 
the problem one might have argued that, once it becomes very strong, the thermoacoustic 
effect may dominate the situation and shift the orientation of the acoustic field to some 
preferable position. This suggests the possibility of the annular prime mover achieving 
onset. Even if the uniform cross section, single stack prime mover does not reach onset, 
other nonuniformities in cross section will have some effect on the spatial distribution of 
the acoustic field. Hence, by placing another constriction in the annulus, the previously 
mentioned pressure and velocity nodes may be displaced from the stack, providing the 
possibility of the annular prime mover reaching onset. These are some of the concepts that 
initiated this work. 

It should be noted that a rigid-rigid prime mover can be considered a limiting case of a 
constricted, annular prime mover. All the previous work on thermoacoustic prime movers 
have dealt with engines that have some sort of well-defined boundary conditions imposed 
upon them. There has been no detailed analysis of prime movers with periodic boundary 
conditions. 



A. BACKGROUND 

Although thermoacoustic phenomena have been observed for a long time, significant 
advances in practical thermoacoustics are relatively recent. The earliest example of a 
thermoacoustic prime mover is the Sondhauss tube (Sondhauss, 1850). Over 100 years 
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ago, it was observed by glassblowers that a hot glass bulb attached to a cool glass tube 
sometimes emitted sound at the tip of the tube. Sondhauss quantitatively investigated the 
relationship between the pitch of the sound emitted and the dimensions of the tube. Lord 
Rayleigh explained the Sondhauss tube quantitatively in 1896, but no complete theoretical 
analysis of these phenomena was made before publication of the series of papers by 
Nikolaus Rott (Rott, 1969, 1980, 1983). Later, Wheatley (1985) and others at Los 
Alamos began the development of practical thermoacoustics devices. 



B. PAST RELEVANT WORKS 

Past works relevant to this research fall under three headings: (1) the thermoacoustic 
prime mover, (2) the acoustic Stirling Engine, and (3) the annular resonator. 

1. The Thermoacoustic Prime Mover 

The potential application of thermoacoustics as a heat-driven sound source has 
motivated theoretical developments (Rott, 1969; Wheatley and Cox, 1988; Swift, 1988). 
Work on prime movers has been focused on either rigid-rigid or rigid-open prime movers. 
Some examples are outlined below. Miglori and Swift (1988) constructed a thermoacoustic 
prime mover that used liquid sodium as its working fluid. A prime mover has also been 
used as a heat-driven sonar projector (Gabrielson, 1991). Swift (1992) built and analyzed 
a 5-inch thermoacoustic prime mover which, at its most powerful operating point, using 
13.8-bar helium, delivered 630 W to an external acoustic load. With minor modifications, 
this device was used by Swift (1994) to research the application of similitude to nonlinear 
thermoacoustics. 
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The Naval Postgraduate School has been in the forefront of much thermoacoustic 
research. Atchley and the author first built and tested a heat driven prime mover (Lin and 
Atchley 1989), which generated a sound, at a temperature difference of 453 °C, with a 
peak acoustic amplitude of 7.9% of atmospheric pressure. Subsequently, other work was 
done on this device (Atchley 1992, 1993; Atchley and Kuo, 1994). In 1993, a 
thermoacoustic prime mover was constructed by Castro and Hofler to investigate the 
performance of heat exchangers and produced peak pressures of up to 20% of the mean 
pressure (Castro and Hofler, 1993). 

DeltaE, a program developed at Los Alamos National Laboratory by Ward and Swift 
(1994, 1996) has been applied to several cases in the design and analysis of thermoacoustic 
prime movers. DeltaE stands for Design Environment for Low Amplitude Thermoacoustic 
Engines. It solves Rott’s wave equation for a geometry defined by the user. For example, 
Swift (1992, 1996) used DeltaE to analyze the performance of a 5-inch thermoacoustic 
engine. Nessler (1994) used DeltaE to compare the performance of the pin stack with a 
conventional stack in a thermoacoustic prime mover (Swift and Keolian 1993). Yang 
(1995) and Meng (1996) also used DeltaE to predict onset of their thermoacoustic prime 
movers. More recently, DeltaE was applied to analyze the efficiency of a thermoacoustic 
prime mover with a pin stack (Gibson, Nessler, and Keolian, 1997). 

2. The Acoustic Stirling Engine 

Ceperley (1979, 1982 and 1985) has discussed acoustic engines using traveling 
waves. As recognized by Ceperley and several other authors (Swift, 1988, 1995; Hofler 
1988; Organ 1992), the phasing between pressure oscillations and velocity within a Stirling 
engine regenerator is the same as those of an acoustic traveling wave. Ceperley proposed 
the replacement of the function of the pistons in a Stirling engine by acoustic processes, to 
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form what he called a traveling-wave heat engine. One of his concepts, a traveling- wave 
heat-driven refrigerator, is shown in Fig. 1.4. In this device, one regenerator and heat 
exchanger set functions as the prime mover, adding acoustic power onto the traveling wave 
as heat flows from a high-temperature heat source to a room temperature heat sink. The 
other set functions as a heat pump, using acoustic power from the traveling wave to 
transport heat from a low temperature to room temperature. His work was partly 
responsible for motivating the author to embark on this project. However, his analysis was 
overly simplistic and he never built a working device. 





Figure 1.4 Ceperley’s design for a traveling-wave heat-driven refrigerator. The 
arrows show the direction of wave propagation. 

3. The Annular Resonator 

Previous work on annular resonators is also relevant to the study of the annular prime 
mover. The eigenmodes for electromagnetic wave propagation in an annular cavity are of 
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fundamental interest in the radio-frequency (RF) heating of tokamak plasmas, besides 
having other applications in the microwave circuit theory (Cap and Deutsch 1978, 1980; 
Janaki and Dasgupta, 1990; Wu, 1992). 

Another category of relevant work is the frequency splitting phenomenon found in 
different fields. As early as 1969, Rudnick et al. made observations of a superfluid helium 
persistent current using the Doppler-shifted splitting of an' azimuthal resonant fourth-sound 
mode of the cylindrical resonator containing liquid helium (Rudnick et al. 1969). 
Heiserman investigated the persistent currents in superleaks in contact with bulk superfluid 
helium using the Doppler shifts of the acoustic modes of an annular resonator partially 
packed with a superleak, and with a simple gyroscopic technique (Heiserman, 1975). In 
plasma physics, in the presence of a finite plasma current, the axisymmetric 
magnetoacoustic wave resonance exhibits a frequency splitting for a finite annular mode 
number between oppositely directed traveling waves (Borg and Wit, 1991). A deformation 
of the cross section of a pinch of the plasma cross section also induces frequency shift and 
splitting of circular modes (Ring, 1988). In the microwave and RF field, frequency 
splitting can also be found in a nonuniform microstrip ring resonator (Wested and 
Andersen, 1972), in an asymmetrically coupled microstrip ring resonators (Al-Charchafchi 
and Dawson, 1990) and in a microstrip rhombic resonator with a step discontinuity in one 
of the arms (Al-Charchafchi and Boulkos, 1990; Al-Charchafchi and Schreck, 1994). 



C. OUTLINE OF THIS WORK 

This work constitutes the first detailed theoretical and experimental investigation of a 
thermoacoustic prime mover with periodic boundary conditions. There are five significant 
aspects to this research: (1) using DeltaE to analyze an annular prime mover, (2) 
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developing an entirely new analysis program using MATLAB, (3) designing, building, and 
experimentally investigating a single stack, annular prime mover, (4) experimentally 
investigating a constricted, single stack annular prime mover, (5) predicting the 
performance of a two stack annular prime mover. 

The major conclusions are: (1) A single stack annular prime mover will not reach 
onset because the eigenmodes of the system do not support thermoacoustic growth. (2) A 
constricted annular prime mover will reach onset because the constriction forces dominating 
boundary conditions that alter the eigenmodes. (3) A two stack prime mover is predicted to 
reach onset because one of the eigenmodes of the symmetric system does support 
thermoacoustic growth. 
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II. THEORY 



This chapter is divided into four parts. The first part is a brief review of the aspects of 
thermoacoustics pertinent to the annular prime mover problem. The review draws heavily 
from Swift’s work (Swift, 1988 and 1997), and is intended to introduce the concepts 
required to understand the techniques used to analyze the annular prime mover. One of the 
central questions to be addressed in this dissertation is the ability of an annular prime mover 
to reach onset. The dependence of quality factor Q of the prime mover on the temperature 

difference AT applied to the stack is a key measure of this ability. In this research, the 

predicted Q will be determined from the complex eigenfrequency of the prime mover. 
Therefore, following the review of thermoacoustics, the complex eigenfrequency will be 
introduced. 

Another important aspect of this dissertation is the numerical analysis of the annular 
prime mover. The techniques used will be discussed in the last part of this chapter. Some 
basic properties of the two-point boundary-value problem are discussed, followed by a 
discussion of the advantages and disadvantages of using DeltaE with the annular prime 
mover geometry. A MATLAB numerical analysis program is then described that is more 
tailored to the annular prime mover than is DeltaE. The validation of this program is 
discussed and the application of the MATLAB program to our annular prime mover is 
presented. 

This chapter concludes with a discussion of the end corrections for constrictions in the 
cross section of annular resonators. 



11 



A. REVIEW OF THERMOACOUSTICS 



Thermoacoustics is essentially the study of the acoustics of a fluid in which there 
exists a temperature gradient. The acoustics of such a fluid is contained in three differential 
equations which relate the complex acoustic pressure and volume velocity, and the 
temperature and enthalpy of the fluid. This method was first described by Wheatley et al. 
(Wheatley et al., 1983). The derivations of these relationships are summarized in this 
section. Once these concepts have been introduced, a description of the basic scheme used 
to analyze the performance of an annular prime mover is given. 

This section is essentially a summary of pertinent aspects of thermoacoustics given in 
the articles by G. W. Swift (Swift, 1988 and 1997). Swift’s quantitative thermoacoustic 
analysis is based on Rott’s wave equation and energy-flux equation, which result from the 
momentum, mass continuity, and energy equations in the acoustic approximation. The 
results are valid for arbitrary phasing between the acoustic pressure and velocity; in other 
words, for both standing-wave and traveling-wave systems. The results presented in this 
section serve as the basis for the numerical analysis used in this dissertation. 

First, the notation that will be used throughout this discussion is established. The 
acoustic wave is considered to propagate in the x direction within a duct of constant cross 
section. The usual complex notation is adopted for time-oscillatory acoustic quantities. 
Also, expansions to first order in the acoustic variables is assumed to suffice. Thus, the 
pressure, velocity and temperature can be written 



p(x) = p m + p l (x) e'“. 


(2.1) 


u(x, y, z) = m, (x, y, z) e“°' , 


(2.2) 



and 
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T(x, ,y, z) = T m (x) + T, (x, y, z) t im , 



(2.3) 



where i = V-T and m, is the x-component of v . The mean values (subscript m) are real, 
but the acoustic variables (subscript 1) are, in general complex, to account for the relative 
phasing of those oscillating quantities. It is assumed that the mean fluid velocity u m = 0. 
All the time dependence appears in the e im term, with (O = 2 n f being the angular 
frequency. For clarity, we consider here only the case of large solid heat capacity, so that 
the temperature of the solid material in the stack is simply T m (x), independent of t, y, and z. 

The conservation of momentum of an incompressible viscous fluid is ( Landau, 1975) 



^ + (v • grad) v = - grad p + ^ V 2 v , (2.4) 

where v is the velocity, p is the density, p is the pressure and /i is the dynamic viscosity of 
the fluid in the resonator. This is called the Navier-Stokes equation. To develop a 
quantitative understanding including viscous effects, we begin by finding the y and z 
dependence of the velocity. Because of the viscous boundary layer in the y and z direction 
all other viscous derivatives can be neglected compared to p9 2 w,/9y 2 and pd 2 u ] /dz 2 . 
Equation (2.4) then reduces to 



icop m u x 



dp, 

dx 



+ fi 



(d 2 u I 

U 2 



d 2 u, 



(2.5) 



Equation (2.5) is an ordinary differential equation for M,(y, z). The boundary 
condition at the solid surface is w, = 0. With this boundary condition Eq. (2.5) can be 
solved to yield 
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( 2 . 6 ) 



«i(*. y, z)= — — [ 1 -h v (y, Z)]^y- , 
cop m dx 



where h v (y, z) depends on the specific geometry under consideration. The volume velocity 
is the spatial average over y and z of Eq. (2.6) and is thus 



Ui(x) = (1-fJ dp ' 



®Pn 



dx 



(2.7) 



where / v is the spatial average of h v (y, z). For the case of parallel-plate geometry of the 
stack used in this work, it can be shown that (Swift 1988; Arnott et al. 1991) 



K (y) = 



cosh[(l + Qy/g v ] 
cosh[(l + Oy 0 /5j’ 



( 2 . 8 ) 



and 



fv = 



tanh[q + oyg v ] 
Kl + i)y 0 /S v ] 



(2.9) 



where S v = ^2jU/p m O) is the fluid’s viscous penetration depth and y 0 is plate half-gap. 

Calculation of the oscillating fluid temperature is accomplished using only the first 
order terms in the general equation of heat transfer (Landau, 1975). The entropy is 
expressed in terms of the acoustic pressure and temperature of the fluid and the resulting 
ordinary differential equation is solved subject to the appropriate boundary conditions. In 
the presence of a temperature gradient along the duct, the solution is (Swift 1988 and 1997) 



T^x, y) = 



1 



P m c p 



[ i AlPi + 



®(l-/v)A gas dx 



dT ~ (i- Vl PS ) Ux , (2.10) 



1 - a 



where <J=c p jl/K is the Prandtl number, and c p is the isobaric heat capacity per unit mass. 
The spatial average in (y, z)-direction of Eq. (2.10) can be written 
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r,(; c) = 



— 1 1 -ap, + 

p.c, 



[ dT m ( 1 f K - g/y 

©(l-/v)V dx l - a 



)U , , 



( 2 . 11 ) 



where / K is the spatial average of / 2 K . The functions /z K and / K are the same as h v and / v , only 
with (5 V replaced by 8 K = ^2k j p m c p (D , the thermal penetration depth, where Kris the thermal 
conductivity of the fluid. The function / K and/ v are very import an tin thermoacoustics and 
describe how viscous and thermal processes affect the oscillatory velocity and temperature 
in the acoustic field. A graph of/ v as a function of y/<5 v is shown in Fig. 2.1. 

The equation of continuity is (Landau, 1975) 



§ + V . (pv) = 0. 

To first order, the x-component of Eq. (2.12) can be reduced to the form 



(2.12) 



ri dx 



(2.13) 



Equation (2.13) can be rewritten in terms of volume velocity 



i co p, + 

Combining Eqs. (2.10) and (2.13) yields 



I d(p„ t/Q 



= 0. 



(2.14) 



dU i 
dx 



id) A 



gas 



Pn,C 



[i+(r-i)/,]p, + 



U-/v) 

(i-/ v )(i-<r)r. * 



(2.15) 



where c is the speed of sound in the gas, and y is the ratio of isobaric to isochoric specific 
heats. 
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We pause to discuss the results derived in the last several paragraphs. Equation (2.7) 
shows that the gas which is much farther than 8 V from the nearest solid surface experiences 
essentially no viscous shear and moves with a velocity that depends only on x. The gas 
that is much closer than S v is nearly at rest. Gas that is between these extremes moves with 
a reduced velocity amplitude and significant phase change. The terms in Eq. (2.11) 
indicate that there are two contributions to the oscillatory temperature. The- first term is 
simply due to the adiabatic compressions and expansions of the fluid. The second term 
comes from the convection of gas parcels along a temperature gradient. The net 
temperature oscillation is just a linear combination of, these two effects. It should be 
pointed out that/ K has the same functional dependence as y/S K as is portrayed in Fig. 2.1 
for/v 

Equation (2.15) has an easy physical interpretation. It indicates that dU\jdx comes 
from a density oscillation in the fluid. This density change can be caused by a pressure 
change or by convection of gas along the temperature gradient. 
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Figure 2.1 The real and imaginary parts and magnitude of f v as functions of the ratio 



y/8 v for the case of parallel plate geometry. 



Next, following the derivation by Swift (1988 and 1997), the time averaged energy 



flux along the stack can be found as 



^ 



Pi 1- 



(1 + <7)(1-/ V *)J 



P m c„ U , 2 



20) A gas (1 - O ’ 2 ) |l - /„ 



dT 

Im (f K + [\as K + A soild * solid) 



dT 

m 

dx 



(2.16) 



where the superscripts * represent the complex conjugate of a complex quantity, and fC soljd 
is the thermal conductivity of the solid. On the basis of Swift’s assumptions, in steady 
state for an engine without lateral heat flows to the surroundings, H 2 along x must be a 
constant throughout the stack. 
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To numerically analyze the performance of thermoacoustic engines, it is recognized 
that Eqs. (2.7), (2.15) and (2.16) comprise a set of three second-order coupled complex 
differential equations in the five variables Re(p,), Im(p,), Re( U i), Im( U i) and T m . These 
equations are the basis of both DeltaE and the MATLAB program presented later in this 
chapter. 

In order to obtain a more useful form of the differential equation for p,(x), Eq. (2.7) 
can be written as 



_L^_ Vi(x) 

dx A gas (l-f v ) 



(2.17) 



In a duct region, the integration is governed by Eqs. (2. 15) and (2. 17) and the temperature 
gradient dTJdx is determined by the measured temperature profile. It is assumed that there 
is no temperature gradient across the heat exchanger. As a result, in a heat exchanger, 
integration is carried out by Eq. (2.17) and a simplification of Eq. (2.15) 



dUx 

dx 



i(0A 



gas 



P m c 






(2.18) 



As for the stack region, integration is based on Eqs. (2.15), (2.17), and (2.19). In 
other words, the time averaged energy flux H 2 is only used in the stack region and is a 
constant throughout the stack region. 

The MATLAB numerical analysis method used in this work starts with specifying 
Re(p,), Im(p,), Re((7i), Im(£/i), T m , Re(/), and Im(/) at a point in the duct. The complex 
pressure and volume velocity just outside the starting end of the ambient heat exchanger can 
be determined by integrating Eqs. (2.15) and (2.17) in the boundary layer approximation 
limit. Next, the complex pressure and volume velocity can be determined just inside the 
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starting end of the ambient heat exchanger by invoking continuity of pressure and volume 
velocity at the end of the ambient heat exchanger. Then, Eqs. (2.17) and (2.18) are used to 
integrate over the ambient heat exchanger. Continuity of pressure and volume velocity is 
used to enter the stack. A value of H l is specified to get dTJdx at the beginning of the 
stack. The solution in integrated to the heat exchanger. Again Eqs. (2.17) and (2.18) are 
integrated over the hot heat exchanger. Finally, the pressure, volume velocity, and the 
temperature at the exit of the hot heat exchanger are matched to those at the starting point. 



B. COMPLEX EIGENFREQUEN C Y 

The transition to onset in a prime mover is conveniently discussed in terms of the 
quality factor Q which can be defined as the ratio of 2k times the energy stored in the 
resonator to the energy dissipated per acoustic cycle. Letting W denote the dissipated 
power, 



Q = (2.19) 

W 

Both W and £ st are second order in the acoustic amplitude. A typical prime mover is 
comprised of five sections: the ambient duct, the ambient heat exchanger, the prime mover 
stack, the hot heat exchanger, and the hot duct. W is then the sum of the power dissipated 
in the five individual sections (Lin, 1989; Atchley, 1992 and 1994) 

W = W a mb+ Wambhx + W pms + W hot hx + Who! • (2.20) 

The subscripts amb, amb hx, pms, hot hx, and hot refer to the ambient duct, the ambient 
heat exchanger, the prime mover stack, the hot heat exchanger, and the hot duct. 
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respectively. The reader should refer to Atchley (Atchley et al., 1992) and Lin (Lin, 1989) 
for a full explanation. 

Instead of Q, we prefer to use the reciprocal of Q, because 1 IQ converges to zero at 
onset rather than diverging to infinity. Also, it is this quantity that is proportional to the 
acoustic power output of the prime mover. It is worth mentioning that 1 IQ is a function of 
the prime mover geometry, the thermophysical properties of the gas, p m , and the 
temperature difference. As the temperature difference along the stack increases from zero, 
the thermal losses in the stack decrease and eventually become negative, representing gain. 
In other words, below onset 1 IQ is positive, at onset II Q is zero and above onset \/Q is 
negative. A complete evolution of Q through onset for a prime mover is discussed by 
Atchley (Atchley, 1994). 

In this work, the Q is determined by the complex eigenfrequency. Complex 
eigenvalues occur when a system has underdamped modes. The excitation of a resonator 
evolves in time as txp(i(O 0 t) where the eigenfrequency is (Swift, 1988; Gamaletsos 1993; 
Amott et al., 1994) 



co 0 = 27tf 0 + i 



: rfo_ 

Q 



( 2 . 21 ) 



where / 0 is the real resonance frequency. The quality factor Q thus can be found by 



Q = 



Re(o? 0 ) 
2Im(a» 0 ) ' 



( 2 . 22 ) 



It should be noted that the complex eigenfrequency approach accounts for power 
generation in the stack and dissipation everywhere in the resonator. If this approach were 
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not used, an artificial acoustic source would have to be introduced into the system to deliver 
power to match boundary conditions (Choe, 1997). 

C. THE ANNULAR GEOMETRY FOR A PRIME MOVER 

As discussed in the previous chapter, the major difference between conventional prime 
mover geometry and that of an annular prime mover is that there are no typical dominating 
boundary conditions in the latter. Instead, an annular prime mover satisfies the periodic 
conditions 



p x (r, 6) = p ] (r, 6 + 2 K), (2.23) 

and 

Ux(r, 6) = Ux(r, 6 + 2n) . (2.24) 



At first glance, this difference does not seem to cause much difficulty in the design and 
analysis of an annular prime mover, until one realizes the fact that almost all the previous 
work on conventional prime movers employ some dominating boundary condition as a 
convenient starting point of their analysis. Elimination of these boundary condition makes 
research on an annular prime mover considerably different from that of a conventional 
prime mover. 

It can be shown that an annular resonator can be modeled as a straight duct having an 
effect length L eff that is related to the inner and outer radii of the annulus and the particular 
mode under consideration (Choe, 1997). Therefore, we will treat the annular prime mover 
as a straight duct having periodic boundary conditions. L eff can be obtained by first writing 
the general solution of the wave equation in cylindrical coordinates which includes a linear 
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combination of the cylindrical Bessel and Neumann functions. Note that the Neumann 
function must be retained because the origin is excluded in the annular resonator. The 
effective circumference of the resonator can be then obtained by applying rigid boundary 
conditions at the inner and outer walls of the annulus. The annulus used in the experiments 
presented later has inner and outer radii of 10.0 and 15.3 cm, respectively. The height of 
the duct is 5.0 cm. The cutoff frequency for cross modes is approximately 6.6 kHz. We 
are interested in the fundamental longitudinal mode (A = L eff ) corresponding to frequencies 
in the 400 ~ 450 Hz range, well below cutoff. 

D. COMPUTER SIMULATION OF THE ANNULAR PRIME MOVER 

This subsection discusses the numerical analysis techniques for the annular prime 
mover. For the purpose of numerical integration in a thermoacoustic engine, Eqs. (2.7), 
(2.15), and (2.16) can be rewritten as a set of five, first order, ordinary differential 
equations in terms of five independent acoustic variables: Re(/?,), Im(p,), Re(£/i), Im(Ui) 



and T m . 



gil.RJJgMe. VxM 

dx L (1 -f v ) 



(2.25) 



ggfei-J i0,p " A ‘-VM 




(2.26) 




= Re 




dx 
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d Im 



d. 




gas_ 

„2 



[i+(r-i)f K ] Pl + 




and 




l"! 2 : ' ' + o/v*) 



(2.29) 





< 



2A ; ( K ( 1 - cj 2 )| 1-/ v | 2 R e (©) 



Our intention is to solve this set of coupled ordinary differential equations for Re(p,), 

Im(/?,), Re( U i ), Im( U 1 ) and T m , subject to the periodic boundary conditions at the ends of 
the prime mover. This is the same basic method employed by DeltaE. 

1. Problem Statement 

When ordinary differential equations are required to satisfy boundary conditions at 
more than one value of the independent variable, the resulting problem is called a two point 
boundary value problem. Determining the deflection of a bar rigidly fixed at both ends is a 
typical example where the conditions specified are the deflections of the elastic curve at the 
supports. Heat-flow problems often fall into this class when temperature or temperature 
gradients are given at two points. The problem encountered here is of this kind, too. 

A major distinction between initial value problems and two point boundary value 
problems, as is pointed out by Press (Press et al., 1992), is that in the former case one is 
able to start an acceptable solution at its beginning (initial values) and just march it along by 
numerical integration to the end (final values). However in the latter case, the boundary 
conditions at the beginning point do not determine an unique solution to start with. For this 
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reason, two point boundary value problems require considerably more effort to solve than 
do initial value problems. Keller (Keller, 1992) has discussed the existence and 
uniqueness theory for two point boundary value problems. 

Two different numerical approaches to the annular prime mover are applied in this 
dissertation: DeltaE and a MATLAB program. 

2. DeltaE 

DeltaE is a computer program, developed by Bill Ward and Greg Swift at Los Almos 
National Laboratory, for modeling and designing thermoacoustic and other one- 
dimensional acoustic apparatus. Basically, DeltaE solves the one-dimensional wave 
equation based on the usual low-amplitude acoustic approximation. It numerically 
integrates the wave equation in a gas or fluid, in a geometry provided by the user as a 
sequence of segments, such as ducts, transducers, compliances, heat exchangers and 
stacks. It uses continuity of oscillating pressure, oscillating volume velocity, and mean 
temperature to pass from the end of one segment to the beginning of the next. It uses the 
appropriate wave equation and temperature equation for each segment. The iteration is 
controlled by global parameters such as frequency and mean pressure, and by local 
parameters such as the geometry of a segment and enthalpy flow. By defining a series of 
inputs, global parameters (the guess vector), the desired outputs, and boundary conditions 
(the target vector), the problem is solved iteratively. Guesses are updated until targets are 
achieved. The number of elements in the target vector must equal the number of elements 
in the guess vector, otherwise the system is over- or under-determined. When the 
computed and specified targets agree to within a specific tolerance, the solution is said to 



24 



converge. All calculations in DeltaE are performed in double precision and the user may set 
the error tolerance for matching targets. 

In general, because a pass of DeltaE’s integration does not reach the end with targeted 
values of all variables, a shooting method is used to adjust chosen initial variables in order 
to hit the target values. The user provides guesses for chosen initial variables. A 
successful convergence of DeltaE depends heavily upon a good choice of guess and target 
members. A good guide to the choice of guess and target variables is discussed by Ward 
and Swift (Ward and Swift, 1996). 

DeltaE is very versatile due to the fact that the elements of the guess vector are not 
limited to the conventional choices consisting of real and imaginary parts of p, and U 1 . 
Any variables that have an effect on the target vector variables can be used as a guess. This 
feature enables DeltaE to compute resonance frequency, temperature, and geometrical 
dimensions in order to satisfy specified boundary conditions. DeltaE is very efficient at 
converging to solutions for some complicated acoustic systems, however, it knows nothing 
about acoustics or physics. For example, it does not recognize that negative frequencies, 
negative pressures or negative lengths are physically meaningless. It simply does the 
integration. Therefore, the reasonableness of the solutions produced will almost always 
depend on the quality of the initial guess vector and choices of the guess-target vector. 

To apply DeltaE to the annular prime mover the boundary conditions are incorporated 
into DeltaE’s target vector. The unknown conditions at the BEGIN segment, which DeltaE 
is supposed to find, are in DeltaE’s guess vector. To implement the periodic boundary 
conditions into DeltaE, a SOFTEND segment is inserted right after the initial BEGIN 
segment in the input file of the DeltaE program. Another SOFTEND segment (or they can 
both be HARDEND segments) is used at the end of the model. Four DIFFTARGETs are 
specified which match the amplitude and phase of p, and U\ at these two SOFTEND 
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segments (Ward, 1997). The SOFTEND segments serve only as a tool to calculate the 
outputs and will be ignored in the target vector. SOFTENDs do not force the impedance to 
be zero. A given temperature across the prime mover stack is achieved by selecting the 
temperature at the hot heat exchanger as a target and the amount of heat input to the ambient 
heat exchanger as a guess vector. We now have five targets and only one guess; four more 
guesses are still required. The other four guesses have been selected to be the amplitude 
and phase of p x and U i at the beginning segment. These guesses and targets (as are shown 
in Appendix B) are summarized in Table 2.1. 



Guess vector 


Target vector 


1. Ip, (begin)l 


1. Ip^begin)! - Ip^end)! < £ 


2. Phase[pi (begin)] 


2. Phasefp^begin)] - phasefp^end)] < £ 


3. 1 U\ (begin)l 


3. \ V\ (begin)l -\U\ (end)l < e 


4. Phase[ U 1 (begin)] 


4. Phase[ U 1 (begin)] - phase[ U 1 (end)] < £ 


5. Heat input at ambient heat exchanger 


5. Hot heat exchanger temperature 



Table 2.1 Summary of the guess and target vectors in the DeltaE input file for the annular 
prime mover, e represents the tolerance for convergence. 



The important parameters for this work are the Q, the resonance frequency and the 
mode shape of the prime mover. To attain these desired quantities, more deliberation has to 
be made in determining the input file. Two methods can be utilized to find the Q with 
DeltaE (Ward, 1997). In the first method, the Q can be obtained indirectly by adding an 
artificial driver into the system, sweeping its frequency, and computing the frequency 
response of the system. The resonance frequency and half power points can be extracted 
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from this analysis. An alternate method is to compute the ratio of the stored energy to the 
dissipated energy. The stored energy is computed by adding an artificial driver, 
determining the pressure and velocity throughout the system and integrating the energy 
density. The dissipated energy is computed from the work done by the driver to excite the 
system. Since DeltaE only computes the pressure and velocity at the exit of a segment, the 
system must be broken into many segments to get a reasonable approximation to the stored 
energy. Thus the frequency response method was chosen to find the Q and resonance 
frequency. A side branch electroacoustic transducer was used as a driver, which has 
frequency independent parameters. A MATHEMATICA program (refer to Appendix A) 
was also developed to compute a least squares fit to the frequency response to find Q and 
resonance frequency. After finding resonance frequency, the system is then driven at 
resonance. To attain a detailed modal shape, the prime mover has to be divided into many 
small segments. The modal shape is a plot of the pressure at the end of each segments of 
the converged solution as a function of segment location. 

DeltaE is a very useful tool to predict how a given thermoacoustic apparatus will 
perform, or for helping the user design an apparatus to achieve a desired performance. 
However, it does have several disadvantages when applied to the annular problem. As was 
mentioned in Chapter I, in the presence of a stack, an annular prime mover exhibits 
frequency splitting. The degree to which a particular mode is excited depends on the 
location of the driver relative to the stack. Because of dissipation, Q for each mode is 
finite, which means that the two modes may overlap in the frequency domain. In terms of 
a driven acoustic system, this manifests itself as a multimode excitation. When only one 
mode is excited, it is straightforward to find Q and resonance frequency by a least squares 
fit to the frequency response curve. However, in the case when the driver excites two 
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mode simultaneously, some other technique, for example pole-zero analysis, must be used 
to determine Q and resonance frequency for each individual mode. 

Another disadvantage comes from trying to map out the mode shape with DeltaE. In 
DeltaE, values of acoustic pressure are only provided at the end of each segment. 
Therefore, to obtain a detailed mode shape, the prime mover must be divided into many 
small segments. This inherent property makes it is laborious and impractical to use DeltaE 
to find the mode shape. 

The other disadvantage is that currently there is no straightforward way to apply an 
externally-imposed temperature gradient on a resonator duct (Ward, 1997). This limitation 
makes it very difficult to model an annular prime mover and incorporate measured 
temperatures along the duct into the numerical model. 

Because of the disadvantages of applying the DeltaE program to our problem a 
decision was made to develop our own program. This program, written in MATLAB, was 
validated by comparing the results for a simplified annular prime mover with DeltaE 
results. 

3. MATLAB Program 
a. Approach 

One important aspect of this dissertation is to develop a numerical analysis 
program that is more tailored to the annular resonator problem than DeltaE. The program 
implements the shooting method because it appears to be a straightforward way to solve 
boundary value problems. This method also provides a systematic approach to taking a set 
of ranging shots that allows us to improve our aim systematically (Keller, 1992; Press et 
al., 1992). The program uses functions that already exist in MATLAB where possible. 
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A complete explanation of the procedure of the shooting method can be found in 
Press (1992) and Gerald (1994). Only a summary is presented here. Most differential 
equations of order higher than first can be reduced to coupled sets of first-order differential 
equations. The two point boundary value problem is equivalent to obtaining the solution to 
a set of N coupled first-order differential equations, satisfying N periodic boundary 
conditions at the starting point x a and at the end point x b . 

We start by writing Eqs. (2.25) through (2.29) in a general form 



with the periodic boundary conditions at the starting point x a and final point x b expressed as 



where, in our problem, N equals five. 

To solve this boundary value problem, an initial value problem is created by 
assuming five initial values. The MATLAB function ode45, which implements fourth and 
fifth order Runge-Kutta formulas is then used to integrate the differential equations from x a 
to x b . The local error term for the fourth-order Runge-Kutta method is O (h 5 ) and the global 
error would be O (h 4 ), where h is the step size of the integration (Celia and Gary, 1992; 
Gerald and Wheatley 1994). Now, at the final point x b , a discrepancy vector F is defined 
with the dimension of N, whose components measure how far the first solution attempt is 
from satisfying the boundary conditions at x b . The problem now becomes one of finding 
the vector F such that 




i = 1,2, , N, 



(2.30) 



?,•(*.. ^ >y*) = y/(* b > yv .y*). * = 12, > n, (2.31) 



F { (y,, y 2 '->y N ) = °> 1 =1 > 2,...,n. 



(2.32) 
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To find the roots of this equation the Newton-Raphson method is employed. This method 
provides a very efficient means of converging to a solution, if a sufficiently good initial 
guess is given (Press, 1992). 

Lety denote the entire vector of values y r In the neighborhood of y , each of the 
functions F x can be expanded in a Taylor series about y 0 as (Press, 1992) 



N 

F i{y) = ^(:vo + $0 = F t{y o) + Zt 1 &c j +0 {^y 2 )- (2.33) 

i=i dy, 

The matrix of partial derivatives, known as the Jacobian matrix, J is defined as 

dF 

J„= -T-- (2.34) 

Sy, 



Since it is difficult to compute the matrix J analytically in this problem, finite differences 
are used to compute this matrix. In the MATLAB program, this is accomplished by 
perturbing each y. individually and finding the value of the vector of AFJAy fpc^). 

Equation (2.34) can be expressed in the matrix form 



F(y + 8y) = F(y) + J*dy +0(5y 2 ) . (2.35) 

To obtain a set of linear equations for the correction vector 8y, terms of order 8y 2 
and higher are neglected. Setting F(y+5y) = 0, gives 

J*Sy = - F(y). (2.36) 

Provided J is nonsingular, <5y is given by 



^ = -r‘F{y)- 



(2.37) 



In practice, J ' 1 is obtained by LU factorization. In MATLAB notation, Sy is given by 
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5y = -J\F(y), 



(2.38) 



where \ is used to denote left division. 

The LU factorization implements Gaussian elimination. Because it involves the 
least amount of arithmetic operations, Gaussian elimination is generally considered to be 
one of the most efficient computation methods in the problem of solving a system of linear 
equations (Leon, 1994; Lindfield and Penny, 1995). It is worth pointing out that using 
J\F(y) instead of inv(J)*F(y) in MATLAB is two or three times faster and produces 
residuals on the order of machine accuracy relative to the magnitude of the data (The Math 
Works, 1994). Once 5y is found, the solution vector is then updated as 

(2-39) 

The process is iterated until the solution converges, i.e. when the value of the normalized 
norm of the discrepancy vector F is smaller than a desired threshold. 

In order to formulate a program to model the annular prime mover, the parameters 
to be updated at each pass of the integrations must be identified. These parameters are 
selected to be: the real and imaginary parts of the volume velocity U\ and the complex 
frequency ft), and the time-averaged energy flux along the stack H 2 . Re(p,) and Im(p,) are 
kept constant at x = 0, which simply normalizes the amplitude of the pressure variation. 
In fact, for given values of Re(p,) and Im(p,), a unique solution can be obtained for a 
certain value of H 2 . The periodic boundary conditions are satisfied by matching the 
acoustic pressure p x and the acoustic volume velocity U i between the starting and final 
points. Also the temperature at the hot heat exchanger is matched to a certain value while 
keeping the temperature of the ambient heat exchanger a constant. 
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b. Description 



The complete MATLAB program is described in Appendix C. This section 
provides a brief explanation of how the program works. First, the coordinate system for 
the prime mover (as is shown in Fig. 2.2) is defined. The hot heat exchanger is at one end 
of the system. The end of the hot heat exchanger is at x = L eff . The choice of the position 
of x = 0 is arbitrary, however by putting the stack and heat exchangers at the end, only 
four distinct regions are required. Had the stack been placed in the middle, there would 
have been five. 
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Figure 2.2 Coordinate system used in the numerical analysis of the annular prime 
mover. 



Initial values of the real and imaginary parts of /?,, a measured temperature profile 
7 m (x) along the duct and ambient heat exchanger, and a temperature for the hot heat 
exchanger, T h are given. The real and imaginary parts of U\, the real and imaginary parts 
of the resonance frequency 0) and the time-averaged energy flux along the stack H 2 are 
guessed. It is assumed that the temperature of the duct and the ambient heat exchanger are 
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held at certain values by external means. A set of reasonable initial guesses, which can be 
estimated from the output of DeltaE, is the key to successful convergence of the program. 
The desired temperature of the hot heat exchanger (T h in the program) is accomplished by 
adjusting the time-averaged energy flux along the stack H 2 . The targets to be matched are 
the real and imaginary parts of p x and U 1 at x = 0 and x = L eff , and the temperature of the 
hot heat exchanger. This matching is accomplished by adjusting the real and imaginary 
parts of U i , the real and imaginary parts the complex angular frequency (0, and the time- 
averaged energy flux along the stack. It is also assumed that there are no temperature 
gradients in the ambient and hot heat exchanger. 

Once the required input parameters are determined. The program solves the 
system of equations for the five unknown variables. If the resulting new p v U i, and T h 
differ significantly from the initial guesses it is necessary to update the guesses and do the 
integration again. During the integration, all the temperature dependent parameters and gas 
parameters are calculated according to the temperature distribution. This whole process is 
repeated to convergence. 

c. Validation 

This program was validated by comparing its solution to that using DeltaE for a 
simple case, in which it is assumed that the temperatures of the entire duct and the ambient 
heat exchanger are held at 293 K by some external means. The desired temperature of the 
hot heat exchanger is achieved by adjusting the quantity H 2 . 

The comparisons of the results from the two programs are shown in Figs. 2.3 
through 2.8. Figure 2.3 shows resonance frequency vs. AT, with AT increasing from 
0 K to 240 K. Figure 2.4 shows \/Q vs. AT over the same temperature difference span 
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as in Fig. 2.3. The overall agreement of resonance frequencies and 1 IQ from the two 
methods is good. In fact, the normalized mean square errors of the resonance frequency 
between two methods are 1.48xl0' 2 % and 6. 10x1 0‘ 3 % for the low and high frequency 
mode, respectively. The normalized mean square error of 1 IQ between two methods are 
0.684% (low frequency mode) and 0.854%(high frequency mode). The value of 1/(2 
calculated from DeltaE is determined from a pole-zero analysis of the frequency response. 
The value of Q so determined is sensitive to the bandwidth and weighting used in the fit. 
The error bar attached to the AT = 200 K data point indicates the uncertainty and accounts 
for most of the disagreement at higher values of AT. One reason for the dependence of Q 
on bandwidth is the multimode excitation. A driver was incorporated in the DeltaE 
program to find the frequency response, from which we obtain the required modal 
information. The driver locations are selected such that at AT = 0 K there is only one 
mode being excited. However, as AT is increased, the mode shapes can be altered and the 
previous driver position may excite both modes simultaneously. This can be seen in Figs. 
2.5 and 2.6 which show the frequency response of the low mode and high mode of the 
prime mover from DeltaE for three values of AT. For higher the AT, the effect of high 
mode excitation is clearly visible as a small bump in the frequency response of the low 
mode. However, no such bump is observed in the frequency response of the high mode. 

Figures 2.7 and 2.8 show the mode shapes for the low frequency mode with 
AT = 0 K and AT = 100 K. Figures 2.9 and 2.10 show the mode shapes for high 
frequency mode with AT = 0 K and AT = 100 K. In these figures, the solid lines 
represent the results of the MATLAB program and the open circles represent the results of 
the DeltaE program. The mode shapes from DeltaE have been normalized to the highest 
value of the pressure amplitude. The amplitude of the mode shape from the MATLAB 
program is determined from a least squares fit to the DeltaE results. 
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It is seen that one of the advantages of the MATLAB program is the smooth 
mode shape curve. The agreement of the mode shapes are good in general. However the 
agreement of the low frequency mode is not as good as for the high frequency mode, 
particularly at high temperature difference. The reasons for this discrepancy at higher 
values of AT are not fully understood. It should be pointed out that the two programs treat 
the heat exchangers differently. The MATLAB program assumes that the mean gas 
temperature at the heat exchanger equals the temperature of the heat exchanger. However, 
DeltaE imposes a thermal impedance so that the gas and heat exchanger temperatures differ. 

Figure 2.8 displays one interesting feature. Both program predict a low standing 
wave ratio. This feature is not present in the high mode. 

Although there remain small disagreements between DeltaE and the MATLAB 
program, the major features are in good agreement. In the chapters to follow, the 
MATLAB program will be used to analyze the annular prime mover under more realistic 
conditions. The major difference is that the measured temperature profiles are incorporated 
in the program. 
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Figure 2.3 Comparison of the calculated resonance frequencies of the annular prime 
mover as determined from the MATLAB program and DeltaE. Symbols are explained 

in the legend. Agreement is excellent at all values of AT. The normalized mean square 
error is 1.48xl0' 2 % for the low mode and 6.10xl0' 3 % for the high mode. 
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Calculated 1/Q VS DeltaT for the annular prime mover 



0.04 



0.035 



+ Low mode: Matlab program 

• Low mode: DeltaE program 

* High Mode: Matlab program 

x High Mode: DeltaE program 



0.03 



+ 



0.025 - 



0.02 



0.015 



0.01 



x X x X 



* * * 



* * « 



0.005 L 



0 



50 



100 150 

Delta T (Kelvin) 



200 



Figure 2.4 Comparison of 1 IQ of the prime mover vs. AT as determined from the 

MATLAB and DeltaE. An error bar is indicated for the low mode at AT = 200 K. 
Because of modal overlap, determination of \IQ in DeltaE is not without uncertainty. 
Agreement is better at low AT. The normalized mean square error is 0.684% for the 
low mode and 0.854% for the high mode. 
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Frequency response of the low mode for the annular prime mover 




Figure 2.5 Frequency response of the prime mover from DeltaE (low mode) at 
AT = 0 K, 100 K and 200 K. Note the small bump at about 436 Hz on 

AT = 100 K and 200K curves which is a result of mode overlap. 
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Frequency response of the high mode for the annular prime mover 




Figure 2.6 Frequency response of the prime mover from DeltaE (high mode) at 
AT = 0 K, 100 K and 200 K. Notice that there is no visible bump from overlap. 



39 



Low mode with DeltaT = 0 Kelvin ; frequency = 419.2+6.038i 




Figure 2.7 Comparison of mode shapes of the prime mover at AT = 0 K for the low 

frequency mode as calculated with the MATLAB program and DeltaE. The beginning 
of the stack region is indicated by the dashed vertical line. 
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Low mode with DeltaT = 100 Kelvin ; frequency = 422.3+6. 192i 




Figure 2.8 Comparison of mode shapes of the prime mover at AT = 100 K for the 

low frequency mode as calculated with the MATLAB program and DeltaE. Although 
both programs predict the same general mode shape, the agreement is not as good at 

AT = 0 K. 
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Normalized pressure 



High mode with dT = 0 Kelvin ; freq = 437.1 +2.01 3i 




Figure 2.9 Comparison of mode shapes of the prime mover at AT = 0 K for the 
frequency mode as calculated with the MATLAB program and DeltaE. 
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High mode with dT = 100 Kelvin ; freq = 437+2. 137i 




high frequency mode as calculated with the MATLAB program and DeltaE. The 
agreement is much better than for the low mode at AT = 100 K. 
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E. END EFFECTS 



The usual boundary conditions applied to a change in cross section are continuity of 
pressure and volume velocity. However, at the discontinuity an additional acoustic 
impedance is introduced owing to the abrupt change of cross sectional area (Miles, 1946 
and 1947; Karal, 1953; Pierce, 1991). As a matter of fact, we have shown (Choe, 1997) 
that it is critical to take into account end corrections at the constriction to get good 
agreement between the measured and predicted resonance frequency and mode shapes for a 
constricted annular resonator. 

Sudden changes in duct cross section produce sudden changes in acoustic pressure 
amplitude. Below cutoff, the effect is equivalent to a lumped-parameter impedance. For 
example, a constriction in a duct gives rise to a concentration of flow through the 
constriction that consequently increases the kinetic energy density. As is discussed by 
Morse and Ingard (1968), if the net increase is concentrated in a region small compared to a 
wavelength, the total increase can be characterized by a lumped acoustic inductance. If the 
constriction also produces an increase in viscous energy loss, this addition can be modeled 
as a lumped acoustic resistance. 

This equivalent lumped impedance can be obtained from various methods. Only two 
methods, the conformal transformation method and the higher order mode method, are 
discussed here. 

After computing the analogous acoustical impedance Z, the boundary condition for the 
pressure becomes. 



p 2 (l) = Pl (l) - Z*U x (l), 



(2.40) 
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where p x is the pressure in the unconstricted region and p 2 is the pressure in the 
constriction. U l is the volume velocity in the unconstricted region and the constriction 
junction is at x = l. The volume velocity boundary condition is unchanged. 



1. Conformal Transformation Method 

The first method used to compute the lumped impedance introduced by a constriction 
is discussed by Morse and Ingard (1968). To obtain the analogous impedance of a 
constriction, they compute the extra kinetic energy and viscous-energy loss produced by 
steady flow through the constriction. To find the steady-state flow, the interior of the half- 
duct is first transformed to the upper half of a complex w plane by the Schwartz-Christoffel 
transformation. Next, depending on the geometry of a specific problem, a further 
transformation is required to go from the w plane to the velocity-potential plane in which 
the lumped impedance can be readily obtained. 

Morse and Ingard applied this technique to the case of a rectangular duct of depth d 
and width b for * < 0 and a for x > 0 (as shown in Fig. 2.1 1). The lumped-circuit 
inductance and resistance corresponding to the sudden increase in the cross section of a 
rectangular duct are 



Notice that both L a and /? a vanish when a = b (i.e., when there is no change in width). 




(2.41) 



and 



pco8 v a-b ( a 2 -b 2 , a + b^ 

— 1 + In 

2 ad b ^ 7iab a-b J 



(2.42) 
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Figure 2.11 Section of a rectangular duct with sudden change of width in the y 
direction. 

2. Higher Order Mode Method 

This method involves including higher order modes near the constriction ends 
(Muehleisen, 1996). In this work, Muehleisen examines the reflection, transmission, 
radiation, and coupling of higher order modes at a discontinuity in finite length rigid wall 
rectangular ducts. Using a method of generalized scattering coefficients, an analytic 
expression for the reflection and transmission of higher modes at size discontinuities is 
developed. 

From the scattering matrix formulation an expression for the lumped acoustic 
impedance can be found (Muehleisen, 1997). This method does not include losses, 
however, and so the resistance obtained by conformal mapping is added to the impedance. 
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III. EXPERIMENTAL APPARATUS AND PROCEDURE 



A. INTRODUCTION 

Descriptions of the experimental apparatus and data acquisition procedure are 
presented in this chapter. The chapter begins with a discussion of the annular resonator. 
This resonator was used in preliminary work to measure the modal properties of a 
constricted annular resonator in the absence of a stack (Choe, 1997). This discussion is 
followed by a detailed description of the annular prime mover. Next, the constricted 
annular prime mover is discussed, followed by a description of how the prime mover is 
instrumented with microphones and thermocouples. The calibration of the microphones is 
discussed in this subsection as well. The chapter concludes with a description of the 
procedure followed during data acquisition. 

B. EXPERIMENTAL APPARATUS 

1. Annular Resonator System 

The annular resonator, shown in Fig. 3.1 (without the top) and Fig. 3.2 (with the top) 
has been used in a preliminary work to investigate the modal properties of a constricted 
annular resonator (Choe, 1997). The inner and outer radius of the resonator is 
10.00 ± 0.01 cm and 15.30 ± 0.01 cm, respectively. The outer side wall is a 1 cm thick 
aluminum ring with 15.30 ± 0.01 cm inner radius. A solid aluminum 10 cm radius 
cylinder that is concentric with the outer ring forms the inner side wall. The height of the 
resonator is 5.00 ± 0.01 cm. The dimensions of the resonator were chosen to give a 
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longitudinal resonance frequency in the 300 to 500 Hz range, which is well below the 
cutoff frequency for the first cross mode, as discussed in Chapter II. 

Since the prime mover will be operated at temperatures over 200 °C, the plexi-glass 
top described in Choe has been replaced by a 0.64 ± 0.01 cm thick aluminum top with a 
radius of 16.01 ± 0.01 cm. The top is bolted with a 10.7 cm long threaded rod at the 
center and three equally-spaced 7.1 cm long clamp bars on the edges. Ten microphones 
are installed in the top to measure the spatial pressure distribution inside the driven 
resonator below onset and observe the behavior of the sound generated above onset. To 
monitor the temperatures in the resonator, 5 type-E thermocouples are installed through the 
top and another 5 type-E thermocouples are installed in the stack/heat exchanger assembly. 
A detailed description of the microphone and thermocouple installation will be presented 
later in this chapter. 

It is very difficult to move the stack and heat exchanger assembly after it has been put 
into the resonator because of the complex wiring and the fragility of the glass components. 
Therefore, three holes (3.50 mm diameter) were drilled in the bottom of the resonator to 
allow the driver location to be moved relative to the stack. Depending on which mode is to 
be excited, the driver was then connected to the desired hole via a 0.4 cm O.D., 3.5 cm 
long plastic tube. The two unused holes were plugged with thumbtacks that had some 
Dow Coming high vacuum grease applied to the tips. 

To ensure tight sealing for a high Q, a thin layer of Dow Coming high vacuum grease 
was applied to the surface of the aluminum ring and the center piece before the top was put 
on and bolted down. In addition to the bolt in the center and the three clamps mentioned 
earlier, 6 extra equally spaced clamp bars are applied around the edge of the top to ensure a 
tight seal. 
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The instrumentation is shown in Fig. 3.3. The acoustic signal is generated by a 40W, 
16 Q. external compression driver. The output of a Hewlett Packard 3562 A dynamic signal 
analyzer was sent through a Techron 5507 power amplifier to the driver. The signal 
analyzer was used to determine the frequency response of the resonator. Typically, the 
resonator was driven with a random noise source of 1 Vrms source level. The pole-zero 
analysis function built into the signal analyzer was then used to determine the Q and the 
resonance frequency. To measure the mode shape, the source of the analyzer was changed 
to a lVrms fixed sinusoidal signal at the resonance frequency found in the previous 
measurement. 

2. Stack Assembly 

The prime mover stack assembly consists of the stack plates and the stack holder. The 
first challenge was to select a suitable material for stack. The stack needed to withstand 
high temperatures and the differential expansion imposed upon it from the high gradients. 
Glass has a low thermal conductivity in comparison to most metals (Johanson Companies, 
1996) and is able to withstand higher temperatures than most plastics and other polymers. 
An AF-45 thin borosilicate glass slide (made by Precision Glass & Optics) was chosen to 
be the material for stack. Some properties of this AF-45 glass are provided in Appendix D. 
AF-45 is a modified borosilicate glass with a high content of BaO and A1 2 0 3 and features a 
low thermal conductivity and a low coefficient of thermal expansion of 45.0 x 10' 7 K'* 
(20-300° C). The glass was accurately cut to a size of: 43.65+/-0.05 mm wide and 
23.90+/-0.10 mm long and is 0. 145+/-0.02 mm thick. 

The stack holder (detailed dimensions of which are given in Appendix E) consists of 
two horizontal bars (5.28 cm long) connected to two vertical bars (4.85 cm long) as shown 
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in Fig. 3.4. Tongue-and-groove joints are machined in the ends of the bars to ensure that 
they form a sturdy frame. Sixty-one 0.508 mm deep, 0.203 mm wide, equally-spaced slits 
were cut in the inside of each of the two horizontal bars. The center-to-center spacing 
between the slits is 0.762 mm. The gaps between the cover glass plates are 0.610 mm. 
After the stack holder was assembled, 61 glass slides were then carefully inserted into the 
slits. The overall porosity for the prime mover stack assembly, including the frame, is 
approximately 0.61. 

3. Ambient Heat Exchanger Assembly 

The ambient heat exchanger is used to remove heat from the stack and to maintain one 
end of the stack close to ambient temperature. The ambient heat exchanger assembly is 
constructed as shown in Fig. 3.5 and the detailed dimensions are given in Appendix E. 
The heat exchanger is designed to ensure that the temperature difference between the center 
and the edge of the heat exchanger is less than one degree C for a 35 W load. The ambient 
heat exchanger holder was made in the same way as the stack holder, except that it uses 
copper instead of aluminum and has only 10 slits. Also, a central bar was incorporated to 
supplement the heat transfer between the heat exchanger and the resonator body. The slits 
are 0.635 mm deep, 0.71 1 mm wide, and equally spaced at 4.064 mm (center to center). 
After the holder was made, twenty Alloy-100, 2.30 mm long, 5.10 mm wide, and 0.68 
mm thick copper fins were cut and soft- soldered into the slits (10 on either side of the 
central bar). The overall porosity for the ambient heat exchanger is 0.62 (including the 
frame). The next step was to soft-solder a 5.256 cm x 4.977 cm rectangular copper woven 
wire cloth (20 x 20 mesh) on one side of the heat exchanger assembly. The function of 
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this wire cloth is to increase the thermal contact between the prime mover stack and the 
ambient heat exchanger. 

The fin spacing is much larger than would be used in a practical prime mover. The 
acoustical properties of an constricted annular resonator are very sensitive to both the 
porosity and flow impedance of the constriction. It was decided to make the flow 
impedance of the heat exchangers small compared to that of the stack. This choice 
necessitated using large fin spacing, which reduced the effectiveness of the heat exchanger. 

4. Hot Heat Exchanger Assembly 

The function of the hot heat exchanger is to supply heat to the stack. It consists of 
nickel-chromium wire woven on a non-conducting frame. Since nickel-chromium wire is 
used to provide heat to the system, the hot heat exchanger has to provide electrical 
insulation to the aluminum resonator. Because of the need for electrical insulation and high 
temperature operation, virgin electrical grade PTFE (Teflon) is chosen for the hot heat 
exchanger frame. The hot heat exchanger assembly is shown in Fig. 3.6. Drawings of 
fabricated parts are provided in Appendix E. The Teflon frame has outer dimensions of 
4.978 cm wide, 5.283 cm high, and 0.813 cm thick. The inside of this frame was milled 
out to a size of 4.470 cm wide and 4.663 cm high. Seventeen holes are drilled onto the 
two vertical sides with an equal spacing of 0.521 cm (8 on one side and 9 on the other 
side). Into each of the holes is inserted a 0-80 stainless steel screw. A 92.1 cm, # 28 gage 
nickel-chromium wire with a total resistance of 12 Q was wound around the 17, 0-80 
screws. To avoid direct contact between the nickel-chromium wire and the stack, stainless 
steel flat washers are used on the screws for spacers. The ends of the nickel-chromium 
wire are wrapped with heat shrinkable flexible PVC tubing to insulate the wire from the 
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resonator. The wire is then fed through Teflon tubing which is placed and epoxied in rhe 
bottom wall of the resonator. 
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Figure 3.1 The annular resonator (without the top). 
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Figure 3.2 The annular resonator (with the top). 
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Figure 3.3 Block Diagram of Instrumentation Setup. 
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Figure 3.4 The prime mover stack assembly. 



56 






Figure 3.5 The ambient heat exchanger assembly. 
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Figure 3.6 The hot heat exchanger assembly. 
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5. Stack/Heat Exchanger Assembly 



The stack assembly, the ambient heat exchanger assembly, and the hot heat exchanger 
assembly were connected with four 5.08 cm long, 00-90 stainless steel threaded rods on 
the two horizontal sides (as shown in Fig. 3.7). Teflon spacers were used to ensure proper 
spacing between the stack and the heat exchangers. The spacing between the glass plate 
and the nickel-chromium wire is approximate 1 mm. Before the stack/heat exchanger 
assembly was placed into the resonator, a piece of thin Teflon tape was wrapped on the 
side of the hot heat exchanger to avoid electrical short circuit. A thin layer of silicone heat 
sink compound was applied to the side of the ambient heat exchanger frame to supplement 
the heat transfer between the ambient heat exchanger and the resonator wall and to seal any 
gaps between the frame and the resonator. 
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Figure 3.7 The stack/heat exchanger assembly. 
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6. Constriction 



Placing another constriction in the resonator provides the possibility for an annular 
prime mover to reach onset. Such constrictions were previously used in related work by 
Choe (Choe, 1997). The constrictions are made of PVC (Ploy vinylchloride) because of its 
rigidity and ease of fabrication. The area ratio is defined by the ratio of open area in the 
constriction (which is the shaded block in Fig. 3.8) to the cross-sectional area of the 
resonator. In Fig. 3.8, S, stands for the open area and S 2 stands for the open cross- 
sectional area of the constriction. The table in Figure 3.8 gives the geometrical parameters 
of the constrictions. Constrictions were used having area ratios of 0.7, 0.3 and 0. 1 with a 
total length of 45°. The 45° length was achieved by stacking two 1 1 .25° constrictions and 
one 22.5° constriction together. Vacuum grease was applied to the joining sides of the 
individual constrictions to ensure a tight fit. 
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Area ratio = S1/S2 
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Figure 3.8(a) Nominal geometry of the constrictions. 
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Figure 3.8(b) Measured geometry of the constrictions. 



7. Microphones 

In this experiment, 8 Panasonic electret WM-60AT microphones and two ENDEVCO 
piezoresistive pressure transducers (Model #8510B-5, Serial # G63T, sensitivity 53.85 
mV/psi, and Model #8530015, Serial # AGLC2, sensitivity 10.89 mV/psi) are used to 
measure the pressure distribution inside the resonator. The two ENDEVCO piezoresistive 
pressure transducers are intended to be used to observe the behavior of the sound generated 
above onset, in case the Pansonic microphones saturate. The Panasonic microphones are 6 

mm in height and 5 mm in diameter, with a sensitivity of -44 dB ± 3 dB re 1 V/Pa. A 

Hewlett Packard 6237 Triple Output Power Supply provided 9V DC to each Panasonic 
microphone via a custom made microphone selector. 

As shown in Fig. 3.9, the 8 Panasonic microphones were spaced by 45 degrees in 
pre-drilled holes (0.605 cm in diameter, 0.50 cm deep). They are glued into the holes with 
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silicone adhesive. At the bottoms of these holes, 0.82 mm diameter holes were drilled to 
expose the microphones to the resonator cavity. The two ENDEVCO piezoresistive 
pressure transducers were flush mounted in the cap through spacers with 0.404 cm in 
diameter threaded holes in the center. 



0 degrees 45 degrees 




Figure 3.9 Installation of 10 microphones on the resonator top. 



Prior to the experiment, the 8 Panasonic microphones were calibrated relative to 
Microphone# 1. The calibration apparatus is shown in Fig. 3.10. A calibrating chamber 
was made from a 19.6 mm long, 21.55 mm ID, and 27.22 mm OD circular PVC pipe. A 
Panasonic micro-speaker (Model #EAS-2P106C) was epoxied into one end of the pipe. 
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The resonator top was removed and the calibration chamber was positioned over each 
microphone. A thin layer of Dow Coming high vacuum grease is applied on the edge of 
the chamber to ensure a tight sealing. The frequency responses were then measured for 
each microphone over a 300 ~ 500 Hz frequency range using a Hewlett Packard 35665A 
Dynamic Signal Analyzer. For the calibration, the signal analyzer was set to swept sine 
mode with a integration time of 20 cycles, a settling time of 20 cycles, and a resolution of 
200 pts/sweep. 
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Figure 3. 10 Calibration of the Panasonic microphones using a small chamber. 

Microphone #1 was selected as the reference for the comparison calibration. The 
frequency response of microphone #1 when covered by the calibration chamber is shown 
in Fig. 3.11. Each microphone’s frequency response was divided by that of microphone 
#1 to create a correction factor. A typical curve for the correction factor is shown in Fig. 
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3.12. After finishing the calibrations, one more measurement was made on Microphone #1 
to ensure the reproducibility of the data. With this particular calibration chamber, a 

difference of less than 1% was achieved. The correction factor for the 8 microphones is 

usually in the range of 0.75 to 1.75. 

When measuring the mode shape, the output from each microphone was divided by 
the appropriate correction factor to obtain the corrected amplitude at a particular location and 
driving frequency. A custom made microphone switch board selects each microphone to 
measure the pressure amplitude at different locations. The microphone output amplitude is 
measured using a Standford Research 530 Lock-in amplifier. 




Figure 3.11 Frequency response of microphone #1 when covered by the small 
calibration chamber. 
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Figure 3.12 Plot of correction factor of microphone #2 referenced to microphone #1. 
The correction factor at 400 Hz is 0.940. 



8. Thermocouple Instrumentation 

Ten Omega type E (Model # 5TC-GC-E-30-72) chromel-constantan 0.254 mm 
diameter glass braid insulated thermocouples are used to monitor the temperatures of the 
stack, the heat exchanger, and inside the resonator. Type E thermocouples were selected 

over other types because they provide the greatest sensitivity (in V/°C) and because of their 

useful temperature range. 

As shown in Fig. 3.13, five thermocouples are used to monitor the temperature of the 
stack and the ambient heat exchanger. Two thermocouples each are placed on the hot and 
ambient sides of the center glass slide, one in the center and the other at the edge. These 
thermocouples are fed through 1 .4 mm diameter holes on the copper woven wire cloth of 
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the ambient heat exchanger. Another thermocouple (thermocouple # 6 in Fig. 3.13) was 
placed on one of the fins of the ambient heat exchanger. The thermocouples were attached 
to the glass plate using JB weld adhesive because its high temperature capability. The size 
of the glue spots is approximate 0.01 mm 2 . These glue spots are approximate 0.5 mm 
from the edge the glass plate. The wires of the five thermocouples were individually 
threaded through 1.4 mm diameter holes in the resonator bottom adjacent to the ambient 
heat exchanger. The other five thermocouples were used to sense temperature inside the 
resonator, as is shown in Fig. 3.14. They were fed through 1.4 mm diameter holes in the 
top and positioned approximately 3 cm inside the resonator. All the holes in the resonator 
were then filled with BIPAX TRA-BOND epoxy to ensure a tight seal. 

The temperature difference across the stack is adjusted by the voltage supplied to the 
hot heat exchanger which is controlled by a Power Design 3650R DC power supply. The 
temperatures inside the resonator and in the stack are measured by a Keithley 740 system 
scanning thermometer. 
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Figure 3.13 Thermocouple placement on the stack and the ambient heat exchanger. 
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Figure 3.14 Thermocouple placement in the resonator top. 

C. EXPERIMENTAL PROCEDURE 

1. Resonator Setup 

The stack/heat exchanger assembly is placed at a fixed location relative to the 
microphones and thermocouples. For the experiment, three driver positions are used: 45°, 
90°, and 180° from the center line of the stack/heat exchanger assembly. Unused driver 
holes are plugged with thumbtacks with some vacuum grease applied to the tips. The 
resonator top is oriented so that 0° is set at the stack assembly center line. The 0° setting is 



69 



performed by rotation of the top by hand. To determine the error in positioning the top, the 
top was positioned five times. It was found that the error by rotation was less than ± 1 .0°. 
So, with the mechanical error ± 0.25° in drilling the microphone holes, the total error of 
microphone positioning is less than ± 1.25°. With the stack/heat exchanger assembly 
positioned, the proper amount of vacuum grease is applied on the open surface of the 
resonator and the top is secured with the clamp bars. 

For the measurements with a constriction in the prime mover, a 45° long constriction 
is inserted at a location that is 90° from the center line of the stack/heat exchanger assembly 
to the center of constriction. 

2. Determining Resonator Characteristics 

After the resonator was set, the Hewlett Packard 3562A Dynamic Signal Analyzer is 
used to characterize the resonator. The swept sine mode of the analyzer provides a sine 
wave signal to the system at a given range of frequencies. The frequency response of the 
system is used to determine the resonance frequency and Q of several modes of the 
resonator. When only one mode is excited, the peak pressure amplitude, frequency, and 
quality factor are same as the those given by pole-zero analysis. When two modes are 
excited simultaneously, pole-zero analysis is used to extract the resonance frequency and 
quality factor. The signal analyzer gives poles in the form a+jb. The resonance is simply 
the value of b and the quality factor can be obtained by b/2a. This method has been 
explained by Choe (1997). These measurements are performed for (1) annular prime 
mover without a constriction, and (2) annular prime mover with three different 
constrictions (45° long, porosities of 0.7, 0.3 and 0.1) placed at 90° relative to the center 
of the stack/heat exchanger assembly. The driver in the first experiment was located at 90° 
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(to excite the low frequency mode) or 180° (to excite the high frequency mode) from the 
center of the stack assembly. In experiment (2), the driver was placed at 45° from the 
center of the stack assembly, between the stack and constriction. 

3. Data Collection 

Data collection is performed in sets of annular prime mover without and with a 
constriction. Before data collection, the output of the Hewlett Packard 6237B Triple output 
power supply is set to DC 9V. The output voltage of the Power Design 3650R power 
supply is also adjusted to a certain value and the temperature at the hot end of the stack is 
allowed to reach steady state (temperature variation < 0.1°C/10 sec). The measurement is 
commenced by recording the temperatures measured by the Keithley 740 System Scanning 
thermometer. The resonance frequency and quality factor is measured by the Hewlett 
Packward 3562A Dynamic Signal Analyzer and the microphone switch board is set to the 
microphone that provides the maximum output. The procedure used to acquire the 
resonance frequency and quality factor was described in subsection 1. After finding the 
resonance frequency, a data run is commenced by changing the source of the Hewlett 
Packward 3562A Dynamic Signal Analyzer to a 1 V pp fixed sine with a driving frequency 
equal to the resonance frequency obtained by previous procedure. The acoustic amplitude 
is measured with each microphone. With the driver running, the switch board selects a 
particular microphone and the value of amplitude at each position is measured with the 
Stanford SR530 Lock-In amplifier. The Lock-In amplifier is set to magnitude-phase 
display mode with a time constant of 1 second. At this time, the temperatures are recorded 
again and same measurement are repeated again. Before the entire set of experiments is 
finished, the final temperatures are recorded. Then, the output voltage of the Power Design 
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3650R power supply is changed to adjust the hot end temperature and the same procedure 
is followed. All data recorded from the 8 microphones are scaled by the correction factor 
obtained from the microphone calibration. The corrected data are then entered in a 
computer program to plot the mode shape. This concludes data acquisition of one mode. 
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IV. RESULTS AND DISCUSSION 



This chapter shows the comparison between theoretical and measured results for two 
series of experiments on the annular prime mover. First, measured temperatures in the 
prime mover are presented. Next, results are presented for resonance frequency, l/<2 vs. 
AT, and mode shape for the unconstricted annular prime mover. Results for the annular 
prime mover with a constriction located 90° from the center of the stack/heat exchanger 
assembly are then discussed. Finally, an analysis of the potential for onset in a two-stack 
annular prime mover is presented. 



A. ANNULAR PRIME MOVER 

The measured resonance frequencies and Q’s are determined as functions of AT from 
pole-zero analysis of the frequency response. The theoretical resonance frequencies and 
quality factors are determined from the complex frequency calculated using the measured 
temperature distribution in the prime mover as an input. Figures 4.1 and 4.2 show the 
measured temperatures as functions of the heater voltage at four locations on the center 
plate of the stack for the low mode and high mode of the prime mover. 

Figure 4.1 shows that the transverse (center-to-edge) temperature difference on the 
ambient side of the stack is at most 5.4 K. However, there is a larger transverse 
temperature profile on the hot side. At the highest heater voltage the hot-side, center-to- 
edge, temperature difference is 14.8 K. The transverse temperature profile leads to some 
uncertainty in the value used for AT. The effects of a transverse AT profile on the 
performance of a stack is difficult to access. Therefore, the transverse profiles will be used 
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to place error bars on the AT values. Figure 4.2 shows the center-to-edge temperature 
difference for the high mode. It is seen that the ambient-side temperatures are essentially 
the same as in Fig. 4.1. However the hot side temperatures show a larger center-to-edge 
difference than in Fig. 4.1. The difference is attributed to the fact that the resonator top had 
to be removed to reposition the driver to switch from exciting the low mode to the high 
mode. It is entirely possible, and observed, that while removing and replacing the top the 
stack/heat exchanger assembly was disturbed. This could have slightly altered the hot heat 
exchanger/stack spacing, resulting in the observed differences. 

In the results to follow, all the low mode data were taken, then the driver was 
repositioned and all the high mode data taken. Therefore, the temperatures displayed in 
Figs. 4. 1 and 4.2 should be representative. Typical measured temperature distributions in 
the prime mover are shown in Figs. 4.3 (the low mode) and 4.4 (the high mode). The 
measured temperature profile in the duct (0 to 0.792 m) and the ambient heat exchanger 
temperature are used as inputs into the MATLAB program. The measured hot end stack 
temperature is a target. 

Figure 4.5 shows the comparison of the measured and computed resonance frequency 
for the low mode and high mode. The driver location is 90° from the center of the 
stack/heat exchanger assembly for the low frequency mode and 180° for the high frequency 
mode. It is seen that there is good agreement between the measured and predicted 
resonance frequencies. The biggest differences between the measured and predicted 
resonance frequencies are approximately 0.86 % and 1.77 % for the low mode and high 
mode, respectively. 

Figure 4.6 shows HQ of the low and high mode for the annular prime mover verses 
the temperature differences across the stack. In general the agreement is good, although 
there is a tendency to under predict the HQ (or over predict Q), indicating the presence of 
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unaccounted losses. For the low mode, l/<2 decreases as the temperature difference is 
increased while for the high mode \IQ increases as the temperature increases. It is evident 
from these results that the annular prime mover will not reach onset under readily attainable 
conditions. As a matter of fact, an extrapolation of the l/<2 curve for the low mode gives an 
onset temperature difference of approximately 1400 K. 

The larger discrepancies of the resonance frequency andT/<2 for the high mode at high 
temperature gradient can be caused by a misjudging of the stack temperature. In the 
original calculations, the temperatures of the hot (T h ) and ambient (T a ) sides of the prime 
mover stack are calculated by a weighted average of the. temperatures at the edge and center 
of either side of the stack 



and 



y _ 2-^hc + T he 

1 u — 
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= 2 r ac + r ae 
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where T hc and T he are the center and edge temperature of the hot side of the stack, 
respectively. T zc and T ae are the center and edge temperature of the ambient side of the stack, 
respectively. If the simple average of T h and T a are used in the computations as opposed to 
the weighted average, the resonance frequencies are almost unchanged (difference < 
0.03%). This negligible change is to be expected since the stack/heat exchanger assembly 
is very short compared to the effective length of the prime mover 
( L stack assembly / L eff = 0.038). As a results, the vertical error bars in Fig 4.5 are smaller 
than the symbols used to plot the data. However, 1 IQ is more sensitive to uncertainties in 
T h and T a than is the resonance frequency since all thermoviscous effects are most 
important in the stack region. The sensitivity to uncertainties in AT is represented by the 
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vertical error bars shown in Fig. 4.6. There is also an uncertainty in the calculation of 
temperature difference across the stack. As shown in Fig. 4.1 and 4.2, the temperatures at 
the center of the hot end and ambient end differ from those at the edge. The largest 
differences at the hot end are 14.8 K and 33.3 K for the low and high mode, respectively 
(heater at 23 V in both cases). The differences of the center and edge temperatures in the 
ambient side are approximately 6 K for both cases. This spread is represented by 
horizontal error bars shown in Figs. 4.5 and 4.6. The other possible source of the 
discrepancies in resonance frequency and 1/(2 at high AT for the high mode is the entire 
temperature distribution of the prime mover. As shown in Figs. 4.3 and 4.4, the major 
change in temperature profile is concentrated at the region adjacent to the hot heat exchanger 
and stack. It is difficult to accurately characterize the entire temperature distribution of the 
prime mover using only 10 thermocouples. This is an important difference between the 
conventional and annular prime movers. In a conventional prime mover, there is typically a 
hot temperature duct and an ambient temperature duct. In an annular prime mover the two 
temperatures must merge along the duct. 

Figures 4.7(a), (b), and (c) show the mode shapes for the low mode at different 
temperature differences. The stack/heat exchanger assembly is located between 
x = 0.768 m and x = 0.792 m, indicated by the dashed line. In this experiment, the 
driver is located 90° from the center of the stack which corresponds to a position of 
0.183 m. The driving frequencies are the resonance frequencies as computed by the 
MATLAB program and the AT’s are calculated using Eqs. 4.1 and 4.2. 

Figure 4.7(a) shows the mode shape for the low mode for AT = 0 K. In this case, the 
stack/heat exchanger assembly acts simply like a constriction and there is a velocity 
antinode near the center of the stack/heat exchanger assembly. The agreement between the 
measured and calculated mode shape is good. 
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Figures 4.7(b) and (c) show the mode shapes for the low mode at temperature 
differences of 80 K and 200 K, respectively. The discrepancy between the measured and 
computed mode shape increases as the temperature difference increases. The discrepancy 
appears as a shift in mode shape. However, in both cases, the measured and predicted 
mode shapes show the same general features. In particular, the standing wave ratio 
decreases with increasing AT. Even with this obvious diminishing in the 'standing wave 
ratio, HQ only changes by 15% from AT = 0 K to AT = 220 K. 

Figure 4.8(a) show the mode shape for the high mode with AT = 0 K. The agreement 
between the measured and calculated mode shape is excellent. As one might expect, there 
is a pressure antinode near the center of the stack/heat exchanger assembly. Figures 4.8(b) 
and (c) show the mode shapes for high mode at temperature difference of 76 K and 192 K, 
respectively. The agreement between measured and calculated mode shapes are better than 
those of the low mode. In contrast to the low mode shapes, the high mode shapes do not 
change appreciably as AT increases. In particular, the standing wave ratio remains high. 
This difference between the two modes may be caused by the discontinuity in the acoustic 
impedance near the end of the stack. For the low mode, there is a velocity antinode near 
the center of the stack which results in a low acoustic impedance. However, the acoustic 
impedance of the stack is much larger for the high mode. As a result, as AT increases, the 
impedance discontinuity should have more effect on the low mode than the high mode. 

The most important result of this section is that the annular prime mover will not reach 
onset under easily attainable conditions The reason is that neither mode shape supports 
thermoacoustic growth. In other words, the low mode has a velocity antinode centered on 
the stack/heat exchanger assembly whereas the high mode has a pressure antinode centered 
on the stack/heat exchanger assembly. It was also found that the standing wave ratio of the 
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low mode diminishes as AT increases. However, this alteration in mode shape does not 
correlate with increased attenuation in the prime mover. 
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Figure 4.1 Temperature of the prime mover stack vs. the heater voltage (low 
mode). 
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Figure 4.2 Temperature of the prime mover stack vs. the heater voltage (high 
mode). 
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Figure 4.3 The measured temperature distribution of the prime mover (low 
mode, heater at 23 V). The stack is located to the right of the dotted line. 
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Figure 4.4 The measured temperature distribution of the prime mover (high 
mode, heater at 23 V). The stack is located to the right of the dotted line. 
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Figure 4.5 Comparison of calculated and measured resonance frequencies vs. AT 
for the annular prime mover. 
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1/Q vs DeltaT for the annular prime mover 




Figure 4.6 Comparison of calculated and measured \IQ vs. AT for the annular 
prime mover. 
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Normalized pressure 



Low mode with dT = 0 Kelvin ; freq = 423.1+5.735i 




Figure 4.7(a) Mode shape for the low mode of the annular prime mover when the 
driver is located 90° from the stack and AT = 0 K. 



85 
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Low mode with dT = 79.63 Kelvin ; freq = 427.4+5.722i 




Figure 4.7(b) Mode shape for the low mode of the annular prime mover when 
the driver is located 90° from the stack and AT = 80 K. 
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Figure 4.7(c) Mode shape for the low mode of the annular prime mover when the 
driver is located 90° from the stack and AT = 200 K. 
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Normalized pressure 



High mode with dT = 0 Kelvin ; freq = 438.3+2.099i 




Figure 4.8(a) Mode shape for the high mode of the annular prime mover when 
the driver is located 180° from the stack and AT = 0 K. 
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High mode with dT = 75.73 Kelvin ; freq = 442.7+2.245i 




Figure 4.8(b) Mode shape for the high mode of the annular prime mover when 
the driver is located 180° from the stack and AT = 76 K. 
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High mode with dT = 192.1 Kelvin ; freq = 453.5+2.962i 




Figure 4.8(c) Mode shape for the high mode of the annular prime mover when 
the driver is located 180° from the stack and AT = 192 K. 
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B. CONSTRICTED ANNULAR PRIME MOVER 



In this section, results for the annular prime mover with a constriction located 90° 
from the center of the stack/heat exchanger assembly are presented. The constriction is 45° 
long and three different porosities 0.7, 0.3, and 0.1 are used. First, the comparisons of the 
measured and calculated resonance frequencies of the constricted prime mover are 
presented. Next, the comparisons of the measured and calculated Q of the constricted 
prime mover are discussed. The end corrections are included in the calculations using the 
two methods discussed in Chapter II. This section concludes with the comparisons of the 
measured and computed mode shapes for the constricted prime mover. 

As mentioned earlier, the measured resonance frequencies and Q’s are determined 
from pole-zero analysis of the frequency response. In the measurement, the driver is 
located 45° from the center of the stack assembly and the constriction. The theoretical 
resonance frequencies and Q’s were calculated for four different cases, represented by A, 
B, C, and D. All the results are listed in Tables 4.1 and 4.2. In these tables, Column A 
represents the calculated results without any end corrections. Columns B and C represent 
the results with end corrections for the constriction using the conformal transformation 
method and higher order mode method, respectively. It is worth pointing out that there are 
no appreciable differences between of the results for HQ calculated from the conformal 
mapping transformation and the higher order mode method when the end corrections in the 
constriction are included. Column D gives the results for the case of applying end 
corrections to both the constriction and the stack using the higher order mode method. The 
end corrections in the prime mover stack are approximated by applying the end corrections 
to one slit and summing the impedance of the slits in parallel. The temperature differences 
were determined by Eqs. (4.1) and (4.2). Table 4.1 lists the calculated and measured 
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resonance frequencies for the low and high modes for the various constrictions. Table 4.2 
lists the calculated and measured Q for the low and high modes for the various 
constrictions. In these tables, Rs represents the area ratio of the constriction, i.e., the 
porosity. 

It is seen that the inclusion of end corrections has significant effects on the predictions 
for the low mode. For example, in the case of Rs = 0.1, the averaged differences 
between the calculated and measured resonance frequencies are 1.94 % (no end 
corrections), 0.67% (end corrections using conformal transformation method), and 0.39% 
(end corrections using higher order mode method). On the contrary, including end 
corrections for the constriction has little effect on the predictions for the high mode. This 
point is evident from Figs. 4.9 and 4.10, where it is seen that inclusion of end corrections 
in the constriction is important for the low mode only. However, inclusions of the end 
corrections in both the constriction and stack has more effect on the high mode than on the 
low mode. 

The point to be taken from these figures is that the program does a good job of 
predicting the resonance frequencies, even in the absence of including end corrections. We 
have confidence in the end corrections applied to the constriction. This confidence is the 
result of preliminary work by Choe (Choe, 1997). Application of end corrections to the 
stack is only intended to be approximate, to indicate that they should have an effect. The 
contribution of this work is to point out the necessity of including end corrections. A 
detailed analysis of stack end corrections is an area for future work. 

Figure 4.11 shows the comparisons of the calculated and measured l/Q of the low 
mode for the unconstricted and constricted prime mover with various porosities. In this set 
of experiments, when the constriction porosity is 0.1, the prime mover reaches onset at a 
temperature difference of 260 K. The horizontal error bar in Fig. 4.11 indicates the 
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uncertainty in determining the onset temperature from the measure temperature distribution 
as discussed previously. In all cases, the calculated results under estimate 1 /Q. Therefore, 
the calculated results show a significantly lower onset value of AT in the case of Rs = 0.1. 
The reason for this discrepancy is not fully understood. Nonetheless, both the measured 
and calculated l/Q indicate the same trends. First, the relative values for l/Q are consistent. 
Also, both the predicted and measured results for the constricted prime mover intersect 
within a narrow range of AT. The predicted values intersect at AT = 100K and the 
measured data intersect at AT = 175 K. The 1 /Q curve for Rs =0.1 begins with the 
highest value of the three constricted cases and ends with the lowest value after the 
intersection at AT =100 K. In contrast, l/Q for Rs = 0.7 begins with the lowest value 
and ends with the largest value at higher AT. l/Q for Rs =0.3 remains in between l/Q 
for the other two porosities. 

The conclusion to be drawn from Fig. 4.11 is that it is predicted that a constricted 
annular prime mover will reach onset, and it does. Although, the agreement between the 
measured and calculated results are only fair, there is very good qualitative agreement. 
This indicates that although there are still some details to be investigated, the general 
behavior is understood. 

Figure 4.12 shows the comparisons between the calculated and measured l/<2 for the 
high mode for the unconstricted and constricted prime mover. It is seen that the 
temperature difference has less effect on l/Q for the high mode than for the low mode. 
Over the temperature span from AT =0 K to AT = 200 K, the change in l/Q is less than 
30 % in all cases. The same comments about the general qualitative agreement apply to the 
high mode results as they do to the low mode results. 

The question to be addressed now is why does a constricted annular prime mover 
reach onset. The answer lies in the effects that a constriction has on mode shapes. 
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Figures 4.13(a) and (b) show the comparison of the measured and computed mode 
shapes of the low mode for Rs = 0. 1 . In this case the constriction is much longer than the 
stack assembly (L constriction /L stack = 3.33) and the porosity of the constriction is much smaller 
than that of the stack assembly. As a result, the mode shape is dominated by the 
constriction. There is a pressure node near the center of the constriction and a pressure 
antinode between the constriction and the stack assembly nearest (but not at) the hot end of 
the stack. End corrections for the constriction are included using the higher order mode 
method. Figure 4.13(a) is for AT = 0 K and Fig. 4.13(b) is for AT =288 K. The 
frequencies are the calculated values. It is observed that there is a good agreement between 
the measured and predicted mode shape at AT =0 K. This mode shape is favorable to 
thermoacoustic growth. The agreement in mode shape worsens at AT =288 K. It should 
be noticed that the imaginary part of the calculated resonance frequency is negative, which 
suggests the onset has been reached. 

Figures 4.14(a) and (b) show the comparison of the measured and computed mode 
shapes of the high mode for Rs = 0.1. The agreement between the measured and 
calculated mode shape is good. As one would expect, there is a pressure antinode near the 
center of the constriction. 

One important common feature of Figs. 4.13 and 4.14 is that the mode shapes are 
essentially determined by the constriction and are independent of the temperature 
difference. This is a key to making the annular prime mover reach onset. 

The main conclusion from this section is that to build a functional annular prime 
mover, the relative locations of the constriction and the stack, and their relative porosities 
are the crucial factors. This is equivalent to incorporating some type of dominating 
boundary conditions into the prime mover to alter the mode shapes. 



94 



Rs=0.1 Low mode 


A 


B 


c 


D 




DeltaT(K) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, mea (Hz) 


0 


300.9 


297.1 


296.2 


295.6 


295.6 


41.1 


302.3 


298.6 


297.6 


296.9 


296.9 


70.8 


303.6 


299.9 


298.9 


298.3 


298.1 


103.9 


305.1 


301.2 


300.3 


299.7 


299.2 


175.4 


309.0 


304.9 


304.2 


303.6 


302.6 


228.2 


312.3 


308 


307.4 


306.8 


305.1 


Rs=0.1 High mode 


A 


B 


C 


D 




DeltaT(K) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, mea (Hz) 


0 


474.2 


474 


474.0 


468.8 


470.5 


41.1 


476.6 


476.5 


476.4 


471.6 


474 


70.8 


478.8 


478.7 


478.6 


473.9 


476.9 


103.9 


481.1 


481 


481.0 


476.4 


479.6 


175.4 


487.7 


487.6 


487.6 


483.2 


486.7 


228.2 


493.3 


493.1 


493.1 


488.9 


491.8 


Rs=0.3 Low mode 


A 


B 


C 


D 




DeltaT(K) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, mea (Hz) 


0 


368.4 


363 


361.7 


361.3 


362 


42.1 


368.6 


363.2 


361.9 


361.6 


361.6 


71.7 


369.3 


363.9 


362.6 


362.3 


361.9 


105.6 


370.4 


364.9 


363.6 


363.4 


362.5 


175.4 


373.4 


367.7 


366.6 


366.3 


364.4 


226.4 


376.8 


370.9 


369.9 


369.6 


366.7 


Rs=0.3 High mode 


A 


B 


C 


D 




DeltaT(K) 


f, cal (Hz) 


f.cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, mea (Hz) 


0 


467.3 


467 


466.9 


461.9 


464.2 


42.1 


468.1 


467.8 


467.7 


463.1 


465.8 


71.7 


469.2 


468.8 


468.7 


464.2 


467.2 


105.6 


470.8 


470.5 


470.4 


466 


469.3 


175.4 


475.3 


474.9 


474.9 


470.7 


473.7 


226.4 


480.0 


479.7 


479.6 


475.5 


478 


Rs=0.7 Low mode 


A 


B 


C 


D 




DeltaT(K) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, mea (Hz) 


0 


419.3 


418.4 


417.9 


417.8 


417.8 


46.8 


419.9 


419 


418.5 


418.4 


418.2 


77.1 


421.9 


421 


420.5 


420.4 


418.8 


112.2 


423.5 


422.5 


422.1 


421.9 


419.9 


179.2 


426.9 


425.9 


425.4 


425.2 


422 


231 


433.9 


429.6 


429.2 


428.8 


424.5 


Rs=0.7 High mode 


A 


B 


c 


D 




DeltaT(K) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, cal (Hz) 


f, mea (Hz) 


0 


438.9 


438.8 


438.7 


434.2 


436.2 


46.8 


442.5 


442.4 


442.3 


438.1 


440.2 


77.1 


445.2 


445.1 


445.0 


440.9 


442.2 


112.2 


448.2 


448.1 


448.0 


444 


444.4 


179.2 


454.2 


454.1 


454.0 


450.2 


448.5 


231 


460.2 


460 


460.0 


456.3 


452 



Table 4.1 Comparisons of the calculated and measured resonance frequencies for the 



constricted annular prime mover. Case A: without end corrections; Cases B, C, D: 
End corrections with different methods. 
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Rs=0.1 Low mode 


A 


B 


C 


D 




DeltaT (K) 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, mea 


0 


50.4 


45.5 


45.4 


38.2 


43.3 


41.1 


68.8 


60.8 


60.9 


48.9 


51.5 


70.8 


92.3 


80.2 


80.2 


61.3 


60.5 


103.9 


148.3 


122.6 


123.4 


84.8 


73.8 


175.4 


-519.6 


-1089 


-1029 


421 


127 


228.2 


-124.3 


-123.3 


-131.9 


-225 


239.5 


Rs=0.1 High mode 


A 


B 


c 


D 




DeltaT (K) 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qhigh, mea 


0 


36.4 


36.4 


36.4 


24.3 


32.6 


41.1 


37.9 


37.9 


37.9 


24.9 


33.7 


70.8 


39.2 


39.1 


39.1 


25.4 


34.2 


103.9 


40.8 


40.8 


40.8 


26 


35.9 


175.4 


45.7 


45.6 


45.7 


27.8 


37.9 


228.2 


51.3 


51.2 


51.2 


29.6 


41.3 


Rs=0.3 Low mode 


A 1 


B 


C 


D 




DeltaT (K) 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, mea 


0 


74.9 


69.6 


69.1 


62.7 


65.1 


42.1 


90.2 


84.8 


83.7 


75.6 


73.8 


71.7 


105.5 


98.6 


98.1 


87.7 


81.1 


105.6 


129.2 


122 


123 


107 


90.9 


175.4 


235.2 


234 


248 


196 


121.7 


226.4 


514.4 


590 


690 


435 


157 


Rs=0.3 High mode 


A 


B 


C 


D 




DeltaT (K) 


Qhigh, cal 


Qhigh, cal 


Qhigh, cal 


Qhigh, cal 


Qhigh, mea 


0 


35.7 


35.7 


35.7 


24 


32.8 


42.1 


36.9 


36.9 


36.9 


24.5 


33.7 


71.7 


37.9 


37.8 


37.8 


24.9 


34.3 


105.6 


39.2 


39.1 


39.1 


25.3 


35.1 


175.4 


42.6 


42.4 


42.4 


26.6 


37.4 


226.4 


45.8 


45.8 


45.8 


27.8 


38.8 


Rs=0.7 Low mode 


A 


B 


c 


D 




DeltaT (K) 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, cal 


Qlow, mea 


0 


98.8 


97.3 


97.2 


94.2 


85.4 


46.8 


108.2 


106.8 


106.6 


103.4 


95 


77.1 


118.4 


116.8 


116.7 


114.4 


98.2 


112.2 


136.7 


134.5 


134.5 


133.6 


103 


179.2 


216.3 


207 


205.8 


210 


114 


231 


312.5 


407.1 


398.2 


444.3 


119.1 


Rs=0.7 High mode 


A 


B 


C 


D 




DeltaT (K) 


Qhigh, cal 


Qhigh, cal 


Qhigh, cal 


Qhigh, cal 


Qhigh, mea 


0 


36.6 


36.6 


36.6 


24.7 


35.2 


46.8 


36.5 


36.5 


36.5 


24.7 


36.4 


77.1 


36.1 


36.1 


36.2 


24.3 


36.3 


112.2 


35.7 


35.7 


35.8 


24 


35.2 


179.2 


34.5 


34.8 


34.9 


23.4 


36.4 


231 


33.7 


34 


34.1 


22.8 


37.2 



Table 4.2 Comparisons of the calculated and measured Q for the constricted 
annular prime mover. Case A: without end corrections; Cases B, C, D: End 
corrections with different methods. 
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Low mode, Resonance frequency vs. DeltaT , Rs = 0.1 




Figure 4.9 Comparisons of the measured and calculated resonance frequencies of 
the low mode for the constricted prime mover with a porosity of 0. 1. Method A 
is the conformal mapping transformation and Method B is the higher order mode. 
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High mode, Resonance frequency vs. DeltaT , Rs = 0.1 




Figure 4.10 Comparisons of the measured and calculated resonance frequencies 
of the high mode for the constricted prime mover with a porosity of 0. 1 . Method 
A is the conformal mapping transformation and Method B is the higher order 
mode. 
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1/Q 



0.04 



1/Q for the low mode of the annular prime mover 
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Figure 4.1 1 Comparisons of the measured and predicted HQ of the low mode for 
the unconstricted and constricted prime mover. Symbol o represents the 

measured HQ of unconstricted prime mover. Symbols *, +, and x represent the 

measured 1 IQ for the constricted prime mover with porosities of 0.1, 0.3, and 
0.7, respectively. The lines represent the predicted values. No end corrections 
are included in the predicted values for the unconstricted prime mover. 
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1/Q for the high mode of the annular prime mover 




Figure 4.12 Comparisons of the measured and predicted HQ of the high mode 
for the unconstricted and constricted prime mover. Symbol o represents the 

measured \IQ of unconstricted prime mover. Symbols *, +, and x represent the 

measured \IQ for the constricted prime mover with porosities of 0.1, 0.3, and 
0.7, respectively. The lines represent the predicted values. No end corrections 
are included in the predicted values for the unconstricted prime mover. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 0 Kelvin ; freq = 296.2+3.259i ; Rs = 0.1 




Figure 4.13(a) Mode shape of the low mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT = 0 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end corrections ; dT = 228.2 Kelvin ; freq = 307.4-1 .165i ; Rs = 0.1 




Figure 4.13(b) Mode shape of the low mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT = 228 K. 

The calculated results are based on the higher order mode method. The frequency 
is the calculated value. 
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Normalized pressure 



High mode with end effect ; dT = 0 Kelvin ; freq = 474+6.514i ; Rs = 0.1 




Figure 4.14(a) Mode shape of the high mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT = 0 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 228.2 Kelvin ; freq = 493.1+4. 817i ; Rs = 0.1 




Figure 4.14(b) Mode shape of the high mode of the constricted prime mover 
(Rs = 0.1) when the driver is located 45° from the stack and AT = 228 K. 

The calculated results are based on the higher order mode method. The frequency 
is the calculated value. 
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C. TWO STACK ANNULAR PRIME MOVER 



The knowledge gained from the last two sections leads to the design of a symmetric 
two stack annular prime mover in which the constriction in the last section is replaced by a 
stack assembly. This proposed design of the annular prime mover, as shown in Fig. 4.15, 
has two identical stack/heat exchanger assemblies. The two stack assemblies are placed 
such that the hot heat exchangers face each other. 

The MATLAB program was modified to predict the performance of this two stack 
annular prime mover, by adding one more guess of the time-averaged energy flux Hi 
along the second stack to match the temperature of the hot heat exchanger. Two different 
cases are studied here. In case 1 the duct between the two heat exchangers is held at 
ambient temperature T a = 293 K while in case 2 it is held at the desired temperature of the 
hot heat exchanger T h . The reason for studying case 1 is that the current apparatus can be 
easily modified to conduct the experiment. However, because of the temperature 
discontinuity near the junction of the duct and hot heat exchanger, it is difficult to model 
this case theoretically and computationally. It is known that the temperature discontinuity 
in the duct will also create a small discontinuity in the work flow (Rott, 1969). On the 
contrary, although the theoretical analysis of case 2 is simpler than case 1 , investigating it 
requires a new apparatus. 

Figure 4.16 shows the calculated resonance frequencies of the two-stack prime mover 
for the two cases. In case 1 , it is seen that the resonance frequencies for the high mode 
increase approximately linearly with temperature difference. However, AT has little effect 
on the resonance frequencies of the low mode in case 1. The resonance frequencies for 
both the low and high mode in case 2 increase approximately linearly with temperature 
difference. Because the duct between the two hot heat exchangers is held at T h , the 
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increase in resonance frequency is bigger than that in case 1. It should be pointed out that 
in case 1 at very low temperature difference (AT < 50 K), the resonance frequencies of 
the low mode and high mode are very nearly equal. The MATLAB program always 
converges to only one mode. For example, at AT = 0 K, the program always converges 
to the high mode, while at AT = 25 K it always converges to the low mode. As AT 
increases, the modes are well separated and no such condition was observed. This 
condition is not present in case 2. It was always possible to distinguish between the low 
mode and high mode. It should also be pointed out that for case 1 at AT = 0 K, the 
resonance frequency for the high mode is 424.6 Hz and is 424.0 Hz for the low mode. 
The results indicate the possibility of mode level repulsion (Pippard, 1985 and 1989; 
Laraza and Denardo 1997). 

Figure 4.17 shows the calculated HQ of the low mode and high mode for the two- 
stack prime mover for the two cases. It is evident that the high mode reaches onset at 
AT = 190 K for case 1. The onset temperature difference for the high mode in case 2 is 
approximately 350 K. An interesting feature different from the HQ for the constricted 
prime mover is the linearly increase and decrease in l/<2 for the low mode and high mode, 
respectively. Similar dependence has been observed for a rigid-rigid prime mover (Lin, 
1989; Atchley, et al., 1992). 

Frequency splitting in the two-stack annular prime mover offers the possibility that the 
harmonics of the fundamental frequency and the overtones of the resonator will not 
coincide, thus detuning the system. This means that if high amplitudes are achieved above 
onset, waveform distortion and shocking should not be a problem (Coppens, 1971; 
Atchley and Hofler, 1990; Atchley et. al. 1993). 

Figures 4.18 to 4.21 show the calculated mode shapes for the low mode and high 
mode for the two cases. As one would expect from symmetry, there is a pressure antinode 
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between the two stacks for the high mode and a pressure node between the two stack 
assembly for the low mode. It is evident that only the high modes in both cases support 
thermoacoustic growth, as is confirmed by Fig. 4.17. 



Hot Heat 
Exchanger (Th) 




Figure 4.15 The two-stack annular thermoacoustic prime mover. 
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Calculated resonance frequency vs. dT for the two stack annular prime mover 
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Figure 4.16 The calculated resonance frequencies for the two stack annular prime 
mover. Case 1: The duct between the two hot heat exchangers is held at T a . 
Case 2: The duct between the two hot heat exchangers is held at T h . 
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1/Q 



Calculated 1/Q vs. dT for the two stack annular prime mover 




Figure 4.17 The calculated 1 IQ for the low mode and high mode for the two 
stack annular prime mover. Case 1: The duct between the two hot heat 

exchangers is held at 7 a . Case 2: The duct between the two hot heat exchangers 
is held at T h . 
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Normalized pressure 



Mode shape for of the two stack prime mover, dT = 200K; freq = 426.1+14.74i 




Figure 4.18 The calculated mode shape of the low mode of the two stack annular 
prime mover at AT = 200 K. Case 1: The duct between the two hot heat 

exchangers is held at T a = 293 K. 
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Normalized pressure 



Mode shape for of the two stack prime mover, dT = 1 99.2K; freq = 433.5-0.61 4i 




Figure 4. 19 The calculated mode shape of the high mode of the two stack annular 
prime mover at AT = 200 K. Case 1: The duct between the two hot heat 
exchangers is held at T a = 293 K. 
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Normalized pressure 



Mode shape for of the two stack prime mover, dT = 200K; freq = 535.9+1 8.2i 




Figure 4.20 The calculated mode shape of the low mode of the two stack annular 
prime mover at AT = 200 K. Case 2: The duct between the two hot heat 
exchangers is held at T h = 493 K. 
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Normalized pressure 



Mode shape for of the two stack prime mover, dT = 199.9K; freq = 540.7+3.941 i 




Figure 4.21 The calculated mode shape of the high mode of the two stack annular 
prime mover at AT = 200 K. Case 2: The duct between the two hot heat 

exchangers is held at T h = 493 K. 
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V. SUMMARY AND CONCLUSION 



The objective of this dissertation is to investigate a thermoacoustic prime mover having 
periodic boundary conditions both experimentally and theoretically. The experimental 
approach has been to design, build, and test an annular thermoacoustic prime mover. The 
resonance frequencies, quality factors, and the mode shapes of the prime mover are 
measured as functions of the temperature difference across the prime mover stack. The 
experimental resonance frequencies and Q’s are determined from pole-zero analysis of the 
measured frequency response. 

An important aspect of this dissertation is the computational analysis of a 
thermoacoustic device with periodic boundary conditions. Two approaches were taken. 
First an existing program, DeltaE, was applied to this novel geometry. It was found that 
DeltaE can be applied to the problem, but has some disadvantages. It is difficult to extract 
Q’s below onset and awkward to extract mode shapes. It is also difficult to include 
temperature gradients in ducts. 

Because of these problems, a new program was developed. It uses the same basic 
approach as DeltaE, but overcomes the disadvantages. The complex eigenfrequency is 
used to extract resonance frequencies and Q’s. 

The results of this dissertation indicate that under all but perhaps extreme conditions an 
annular prime mover will not reach onset. The reason is well understood. The eigenmodes 
of a single-stack annular prime mover do not support thermoacoustic growth. This finding 
lead to the investigation of a constricted annular prime mover. The idea is that the 
constriction could impose dominating boundary conditions on the prime mover thus 
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altering the eigenmodes. It was shown that decreasing the porosity of the constriction 
forced the prime mover into onset. The porosity and location of the constriction is used to 
adjust the eigenmodes of the system so that thermoacoustic growth is supported. 

Another way to produce eigenmodes that support thermoacoustic growth is to 
introduce a symmetry into the system which shifts the nodes and antinodes of the modes 
off the stacks, such that for one of the modes a pressure antinode is near the hot end of the 
stacks. It is predicted that a two-stack annular prime mover will reach onset. This design 
is an area for future research. 

Another area for future study is end correction for stacks. It was shown in this 
dissertation that inclusion of end corrections does affect the predicted resonance frequencies 
and Q’s. Determining the proper end corrections for realistic stack geometries is likely to 
be a nontrivial project. 
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APPENDIX A.l THE DELTAE INPUT FILE FOR THE LOW FREQUENCY 
MODE OF THE ANNULAR THERMOACOUSTIC PRIME MOVER 



TITLE Annular resonator (file : to2dTflowC.in) 

!Created@23: 18:35 02-MAY-97 with DeltaE Vers. 3.0bl for the Power Macintosh 
! 0 



BEGIN Initial 












1.0130E+05 


a Mean P 


Pa 




16.63 


A lpl@0 G(0d) 


P 


427.0 


b Freq. 


Hz 




-43.90 


B Ph(p)0 G( Oe) 


P 


293.0 


c T-beg 


K 




6.7062E-04 


CIUI@0 G( Of) 


P 


16.63 


d lpl@0 


Pa 


G 


38.50 


D Ph(U)0 G( Og) 


P 


-43.90 


e Ph(p)0 


deg 


G 


-1.8309E-03 


E Heatln G(27e) 


P 


6.7062E-04 


f IUI@0 


m A 3/s 


G 








38.50 


g Ph(U)0 


deg 


G 








air Gas 


type 












ideal 

i 


Solid type 


1 










SOFTEND 


1 


1 










0.0000 


a Re(Z) 


(t) 




16.63 


A Ipl Pa 




0.0000 


b Im(Z) 


(t) 




-43.90 


B Ph(p) deg 












6.7062E-04 


C IUI m A 3/s 












38.50 


D Ph(U) deg 












7.3673E-04 


E Hdot W 












7.3673E-04 


F Work W 












2.0847E-02 


G Re(Z) 




sameas 0 


Gas type 






-0.1564 


H Im(Z) 




ideal 

| 


Solid type 


o 




293.0 


I T K 




ISODUCT 


Duct 












2.6300E-03 


a Area 


m A 2 




42.05 


A Ipl Pa 




0.2050 


b Perim 


m 




-48.73 


B Ph(p) deg 




3.1788E-02 


c Length 


m 




6.2404E-04 


C IUI m A 3/s 












38.20 


D Ph(U) deg 




sameas 0 


Gas type 






7.0088E-04 


E Hdot W 
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ideal 

i 


Solid type 


7 


7.0088E-04 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


64.95 


A Ipl 


Pa 


0.2050 


b Perim 


m 


-49.96 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


5.3889E-04 


C IUI 


m A 3/s 








37.86 


D Ph(U) 


deg 


sameas 0 


Gas type 




6.6858E-04 


EHdot 


W 


ideal 

i 


Solid type 


A 


6.685 8E-04 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


83.85 


A Ipl 


Pa 


0.2050 


b Perim 


m 


-50.55 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


4.2045E-04 


C IUI 


m A 3/s 








37.36 


D Ph(U) 


deg 


sameas 0 


Gas type 




6.4101E-04 


EHdot 


W 


ideal 


Solid type 


c 


6.4101E-04 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


97.57 


A Ipl 


Pa 


0.2050 


b Perim 


m 


-50.93 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


2.7606E-04 


C IUI 


m A 3/s 








36.44 


D Ph(U) 


deg 


sameas 0 


Gas type 




6.1820E-04 


EHdot 


W 


ideal 

f 


Solid type 




6. 1 820E-04 


F Work 


W 


ISODUCT 


Duct 


— o — 








2.6300E-03 


a Area 


m A 2 


105.3 


A Ipl 


Pa 


0.2050 


b Perim 


m 


-51.21 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


1.1484E-04 


C IUI 


m A 3/s 








33.10 


D Ph(U) 


deg 


sameas 0 Gas type 




5.9903E-04 


EHdot 


W 


ideal Solid type 

I 


n 


5.9903E-04 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


106.4 


A Ipl 


Pa 
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0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 




ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 




ISODUCT 


Duct 7 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


i 


ISODUCT 


Duct 


8 


2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 
ideal 

l 

ISODUCT 


Gas type 
Solid type 


i 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 





-51.45 


B Ph(p) 


deg 


5.5371E-05 


C IUI 


m A 3/s 


-130.1 


D Ph(U) 


deg 


5.8148E-04 


EHdot 


W 


5.8148E-04 


F Work 


W 



101.0 


A Ipl 


Pa 


-51.67 


B Ph(p) 


deg 


2.1977E-04 


C IUI 


m A 3/s 


-138.8 


D Ph(U) 


deg 


5.6316E-04 


EHdot 


W 


5.6316E-04 


F Work 


W 



89.39 


A Ipl 


Pa 


-51.92 


B Ph(p) 


deg 


3.7131E-04 


C IUI 


m A 3/s 


-140.0 


D Ph(U) 


deg 


5.4183E-04 


EHdot 


W 


5.4183E-04 


F Work 


W 



72.22 


A Ipl 


Pa 


-52.23 


B Ph(p) 


deg 


4.9997E-04 


C IUI 


m A 3/s 


-140.6 


D Ph(U) 


deg 


5.1603E-04 


EHdot 


W 


5.1603E-04 


F Work 


W 



50.59 


A Ipl 


Pa 


-52.75 


B Ph(p) 


deg 


5.9773E-04 


C IUI 


m A 3/s 


-140.9 


D Ph(U) 


deg 


4.8535E-04 


EHdot 


W 
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ideal 


Solid type 


4.8535E-04 


F Work 


W 


I 

ISODUCT 

2.6300E-03 


12 — 

Duct 

a Area m A 2 


25.84 


A Ipl 


Pa 


0.2050 


b Perim m 


-54.17 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


6.5855E-04 


C IUI 


m A 3/s 


sameas 0 


Gas type 


-141.1 

4.5059E-04 


D Ph(U) 
EHdot 


deg 

W 


ideal 


Solid type 


4.5059E-04 


F Work 


W 


1 

ISODUCT 

2.6300E-03 


13 — 

Duct 

a Area m A 2 


1.340 


A Ipl 


Pa 


0.2050 


b Perim m 


-165.9 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


6.7865E-04 


CIUI 


m A 3/s 


sameas 0 


Gas type 


-141.3 

4.1353E-04 


D Ph(U) 
EHdot 


deg 

W 


ideal 


Solid type 


4.1353E-04 


F Work 


W 


i 

ISODUCT 

2.6300E-03 


14 - 

Duct 

a Area m A 2 


26.91 


A Ipl 


Pa 


0.2050 


b Perim m 


131.0 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


6.5679E-04 


C IUI 


m A 3/s 


sameas 0 


Gas type 


-141.5 

3.7652E-04 


DPh(U) 

EHdot 


deg 

W 


ideal 


Solid type 


3.7652E-04 


F Work 


W 


i 

ISODUCT 

2.6300E-03 


15 ... 

Duct 

a Area m A 2 


51.56 


A Ipl 


Pa 


0.2050 


b Perim m 


129.7 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


5.9432E-04 


CIUI 


m A 3/s 


sameas 0 


Gas type 


-141.6 

3.4190E-04 


D Ph(U) 
EHdot 


deg 

W 


ideal 


Solid type 


3.4190E-04 


F Work 


W 


i 

ISODUCT 

2.6300E-03 


16 - 

Duct 

a Area m A 2 


73.03 


A Ipl 


Pa 
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0.2050 

3.1788E-02 



b Perim 
c Length 



m 



m 



sameas 0 


Gas type 




ideal 

! 


Solid type 


.... 17 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


— 18 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

t 


Solid type 


.... 19 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


— 20 


VDUCER 


Driver 




1.0000E-09 


a Re(Ze) 


ohms 


0.0000 


b Im(Ze) 





1 .0000E+04 c Re(T 1 ) V-s/m A 3 
0.0000 d Im(Tl) V-s/m A 3 

-1.0000E+04 e Re(T2) Pa/A 



129.2 


B Ph(p) 


deg 


4.951 IE-04 


C IUI 


m A 3/s 


-141.8 


D Ph(U) 


deg 


3.1142E-04 


EHdot 


W 


3.1 142E-04 


F Work 


W 



89.99 


A Ipl 


Pa 


128.9 


B Ph(p) 


deg 


3.6528E-04 


C IUI 


m A 3/s 


-142.1 


D Ph(U) 


deg 


2.8583E-04 


EHdot 


W 


2.8583E-04 


F Work 


W 



101.4 


A Ipl 


Pa 


128.8 


B Ph(p) 


deg 


2.1289E-04 


C IUI 


m A 3/s 


-142.6 


D Ph(U) 


deg 


2.6468E-04 


EHdot 


W 


2.6468E-04 


F Work 


W 



106.5 


A Ipl 


Pa 


128.7 


B Ph(p) 


deg 


4.7479E-05 


C IUI 


m A 3/s 


-146.9 


D Ph(U) 


deg 


2.4644E-04 


EHdot 


W 


2.4644E-04 


F Work 


W 



106.5 


A Ipl 


Pa 


128.7 


B Ph(p) 


deg 


1.4217E-04 


CIUI 


m A 3/s 


-169.5 


D Ph(U) 


deg 


3.5727E-03 


EHdot 


W 
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0.0000 


f Im(T2) 


Pa/A 


1.0000E-09 


g Re(Zm) Pa-s/m A 3 


1.0000E-09 


h Im(Zm) 


Pa-s/m A 3 


1.000 


i AplVol 


V 


sameas 0 


Gas type 




ideal 

i 


Solid type 


— 21 - 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


22 - 


SODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


— 23 - 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


— 24 - 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 



3.5727E-03 


F Work 


W 


3.3263E-03 


G Workln 


W 


1.000 


H Volts 


V 


1.0651E-02 


I Amps 


V 


-51.35 


J Ph(Ze) 


deg 


1.0000E-04 


K IUxI m A 3/s 


-180.0 


L Ph(-Ux 


deg 



108.1 


A Ipl 


Pa 


127.3 


B Ph(p) 


deg 


7.9112E-05 


CIUI 


m A 3/s 


93.49 


D Ph(U) 


deg 


3.5543E-03 


EHdot 


W 


3.5543E-03 


F Work 


W 



103.1 


A Ipl 


Pa 


125.9 


B Ph(p) 


deg 


2.2143E-04 


CIUI 


m A 3/s 


53.91 


D Ph(U) 


deg 


3.5352E-03 


EHdot 


W 


3.5352E-03 


F Work 


W 



91.74 


A Ipl 


Pa 


124.2 


B Ph(p) 


deg 


3.7161E-04 


C IUI 


m A 3/s 


46.13 


D Ph(U) 


deg 


3.5132E-03 


EHdot 


W 


3.5132E-03 


F Work 


W 



74.83 


A Ipl 


Pa 


122.0 


B Ph(p) 


deg 


5.0174E-04 


CIUI 


m A 3/s 
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sameas 0 


Gas type 


42.72 

3.4867E-03 


D Ph(U) 
EHdot 


deg 

W 


ideal 


Solid type 


3.4867E-03 


F Work 


W 


! 

ISODUCT 

2.6300E-03 


— 25 — 

Duct 

a Area m A 2 


53.48 


A Ipl 


Pa 


0.2050 


b Perim m 


118.2 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


6.0191E-04 


CIUI 


m A 3/s 


sameas 0 


Gas type 


40.63 

3.4553E-03 


D Ph(U) 
EHdot 


deg 

W 


ideal 


Solid type 


3.4553E-03 


F Work 


W 


I 

ISODUCT 

2.6300E-03 


26 — 

Duct 

a Area m A 2 


29.40 


A Ipl 


Pa 


0.2050 


b Perim m 


108.6 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


6.6543E-04 


C IUI 


m A 3/s 


sameas 0 


Gas type 


39.08 

3.4195E-03 


D Ph(U) 
EHdot 


deg 

W 


ideal 


Solid type 


3.4195E-03 


F Work 


W 


i 

HXFRST 
sameas 3a 


27 — 

Ambient HX 
a Area m A 2 


23.02 


A Ipl 


Pa 


0.6290 


b GasA/A 


103.2 


B Ph(p) 


deg 


5.0800E-03 


c Length m 


6.6930E-04 


CIUI 


m A 3/s 


1.7018E-03 


d yO m 


38.94 


D Ph(U) 


deg 


-1.8309E-03 


e Heatln W G 


1.5886E-03 


EHdot 


W 


293.0 


f Est-T K (t) 


3.3455E-03 


F Work 


W 


sameas 0 


Gas type 


-1.8309E-03 


GHeat 


W 


copper 


Solid type 


293.0 


H MetalT K 


i 

STKSLAB 

2.6300E-03 


28 — 

Prime Mover Stack 
a Area m A 2 


16.26 


A Ipl 


Pa 


0.6220 


b GasA/A 


-43.67 


B Ph(p) 


deg 


2.4000E-02 


c Length m 


6.7080E-04 


C IUI 


m A 3/s 


3.0480E-04 


d yO m 


38.51 


D Ph(U) 


deg 


7.5000E-05 


e Lplate m 


1.5886E-03 


EHdot 


W 
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7.4220E-04 


F Work 


W 










293.0 


G T-beg 


K 


sameas 0 


Gas type 






293.0 


H T-end 


K 


kapton 

! 


Solid type 




29 


-2.603 3E-03 


I StkWrk 


W 


HXLAST Cold HX 












sameas 3a 


a Area 


m A 2 


16.63 


A Ipl Pa 


0.7050 


b GasA/A 






-43.90 


B Ph(p) 


deg 


3.0480E-04 


c Length 


m 




6.7062E-04 


C IUI m A 3/s 


1.2540E-03 


d yO 


m 




38.51 


D Ph(U) 


deg 


0.0000 


e Heatln 


w 


(t) 


7.3668E-04 


EHdot 


W 


293.0 


f Est-T 


K 


=29H? 


7.3668E-04 


F Work 


W 


sameas 0 


Gas type 






-8.5195E-04 


G Heat 


W 


copper 

i 


Solid type 




30 — - 


293.0 


H MetalT 


K 


SOFTEND 














0.0000 


a Re(Z) 




(t) 


16.63 


A Ipl Pa 


0.0000 


b Im(Z) 




(t) 


-43.90 


B Ph(p) 


deg 










6.7062E-04 


C IUI m A 3/s 










38.51 


D Ph(U) 


deg 










7.3668E-04 


EHdot 


W 










7.3668E-04 


F Work 


W 










2.0846E-02 


G Re(Z) 




sameas 0 


Gas type 






-0.1564 


H Im(Z) 




ideal 


Solid type 






293.0 


I T K 


i 






31 — - 








DIFFTAR 


Ipl mismatch 










0.0000 


a TargDi 




=31 A? 


8.0065E-05 


A D1-D2 




1A b DIAddr 














30 A c D2Addr 












I 






32 — - 








DIFFTAR 


Ph(p) mismatch 








0.0000 


a TargDi 




=32A? 3.3780E-04 


A D1-D2 




IB b DIAddr 














3 OB c D2Addr 













33 
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DIFFTAR IUI mismatch 

0.0000 a TargDi =33A? 1.3884E-09 A D1-D2 

1C b DIAddr 
30C c D2Addr 

t 34 

DIFFTAR Ph(U) mismatch 

0.0000 a TargDi =34A? -1.0691E-04 A D1-D2 

ID b DIAddr 
30D c D2Addr 

! The restart information below was generated by a previous run 
! You may wish to delete this information before starting a run 
! where you will (interactively) specify a different iteration 
! mode. Edit this table only if you really know your model! 
INVARS 504050607 27 5 
TARGS 5 29 6 31 1 32 1 33 1 34 1 
SPECIALS 0 

PLTVAR 702-1 1010203040571 20 1 
! Plot start, end, and step values. May be edited if you wish. 

! Outer Loop: I Inner Loop . 

410.0 427.0 0.2000 1.000 1.000 1.000 
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APPENDIX A.2 THE DELTAE INPUT FILE FOR THE HIGH FREQUENCY MODE 
OF THE ANNULAR THERMOACOUSTIC PRIME MOVER 



TITLE Annular resonator (file :: to2dTfhighC.in) 

!Created@ 12:00:08 02-MAY-97 with DeltaE Vers. 3.0bl for the Power Macintosh 

! 0 - 

BEGIN Initial 



1.0130E+05 


a Mean P 


Pa 




0.0000 


A lpl@0 G(0d) 


P 


430.5 


b Freq. 


Hz 




0.0000 


B Ph(p)0 G( Oe) 


P 


293.0 


c T-beg 


K 




0.0000 


C IUI@0 G( oo 


P 


500.0 


d lpl@0 


Pa 


G 


0.0000 


D Ph(U)0 G( Og) 


P 


0.0000 


e Ph(p)0 


deg 


G 


0.0000 


E Heatln G(27e) 


P 


2.9000E-04 


f IUI@0 


m A 3/s 


G 








0.0000 


g Ph(U)0 


deg 


G 









air Gas type 

ideal Solid type 

! 1 

SOFTEND 1 



0.0000 a Re(Z) 




(o 


0.0000 


A Ipl 


Pa 


0.0000 b Im(Z) 




(o 


0.0000 


B Ph(p) 


deg 








0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 








0.0000 


EHdot 


W 








0.0000 


F Work 


W 








0.0000 


G Re(Z) 




sameas 0 Gas type 






0.0000 


H Im(Z) 




ideal Solid type 

t 




2 — 


0.0000 


I T 


K 


ISODUCT Duct 












2.6300E-03 a Area 


m A 2 




0.0000 


A Ipl 


Pa 


0.2050 b Perim 


m 




0.0000 


B Ph(p) 


deg 


3.1788E-02 c Length 


m 




0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 Gas type 






0.0000 


EHdot 


W 


ideal Solid type 






0.0000 


F Work 


W 
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3 



ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 


Solid type 


A 


0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 


Solid type 




0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


CIUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 

t 


Solid type 




0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 


Solid type 


7 


0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 
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3.1788E-02 



c Length m 



sameas 0 


Gas type 




ideal 

i 


Solid type 


< 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

i 


Solid type 


i 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 

| 


Solid type 


i 


ISODUCT 


Duct 




2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 
i. . 


Solid type 


— 1 


ISODUCT 


Duct 


2.6300E-03 


a Area 


m A 2 


0.2050 


b Perim 


m 


3.1788E-02 


c Length 


m 


sameas 0 


Gas type 




ideal 


Solid type 





0.0000 


CIUI 


m A 3/s 


0.0000 


D Ph(U) 


deg 


0.0000 


EHdot 


W 


0.0000 


F Work 


W 



0.0000 


A Ipl 


Pa 


0.0000 


B Ph(p) 


deg 


0.0000 


CIUI 


m A 3/s 


0.0000 


D Ph(U) 


deg 


0.0000 


EHdot 


W 


0.0000 


F Work 


W 



0.0000 


A Ipl 


Pa 


0.0000 


B Ph(p) 


deg 


0.0000 


CIUI 


m A 3/s 


0.0000 


D Ph(U) 


deg 


0.0000 


EHdot 


W 


0.0000 


F Work 


W 



0.0000 


A Ipl 


Pa 


0.0000 


B Ph(p) 


deg 


0.0000 


C IUI 


m A 3/s 


0.0000 


D Ph(U) 


deg 


0.0000 


EHdot 


W 


0.0000 


F Work 


W 



0.0000 


A Ipl 


Pa 


0.0000 


B Ph(p) 


deg 


0.0000 


C IUI 


m A 3/s 


0.0000 


D Ph(U) 


deg 


0.0000 


EHdot 


W 


0.0000 


F Work 


W 
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I 



12 



ISODUCT 


Duct 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


CIUI 


m A 3/s 






0.0000 


DPh(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


! 


13 .... 








ISODUCT 


Duct 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


CIUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


! 


14 — 








VDUCER 


Driver 








1.0000E-09 


a Re(Ze) ohms 


0.0000 


A Ipl 


Pa 


0.0000 


b Im(Ze) 


0.0000 


B Ph(p) 


deg 


1.0000E+04 


c Re(Tl) V-s/m A 3 


0.0000 


C IUI 


m A 3/s 


0.0000 


d Im(Tl) V-s/m A 3 


0.0000 


D Ph(U) 


deg 


-1.0000E+04 e Re(T2) Pa/A 


0.0000 


EHdot 


W 


0.0000 


f Im(T2) Pa/A 


0.0000 


F Work 


W 


1.0000E-09 


g Re(Zm) Pa-s/m A 3 


0.0000 


Ci Workln W 


1.0000E-09 


h Im(Zm) Pa-s/m A 3 


0.0000 


H Volts 


V 


1.000 


i AplVol V 


0.0000 


I Amps 


V 






0.0000 


J Ph(Ze) 


deg 


sameas 0 


Gas type 


0.0000 


K lUxI 


m A 3/s 


ideal 


Solid type 


0.0000 


L Ph(-Ux deg 


t 


15 — 








ISODUCT 


Constriction 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


C IUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 
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sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


1 


- 16 — 








ISODUCT 


Duct 7 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


C IUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


1 


17 - 








ISODUCT 


Duct 8 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


CIUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


i 


18 - 








ISODUCT 


Duct 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


C IUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


i 


19 ... 








ISODUCT 


Duct 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


C IUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


i 


20 - 








ISODUCT 


Constriction 
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2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


CIUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 

i 


Solid type 


■— 21 — 


0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 

i 


Solid type 


.— 22 - 


0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 

t 


Solid type 


23 - 


0.0000 


F Work 


W 


ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 




0.0000 


EHdot 


W 


ideal 


Solid type 




0.0000 


F Work 


W 


! 




— 24 - 








ISODUCT 


Duct 










2.6300E-03 


a Area 


m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim 


m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length 


m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) 


deg 
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sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 


1 


----- 25 - 








ISODUCT 


Duct 








2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 


0.2050 


b Perim m 


0.0000 


B Ph(p) 


deg 


3.1788E-02 


c Length m 


0.0000 


C IUI 


m A 3/s 






0.0000 


D Ph(U) 


deg 


sameas 0 


Gas type 


0.0000 


EHdot 


W 


ideal 


Solid type 


0.0000 


F Work 


W 
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ISODUCT Duct 



2.6300E-03 


a Area m A 2 


0.0000 


A Ipl 


Pa 




0.2050 


b Perim m 


0.0000 


B Ph(p) deg 




3.1788E-02 


c Length m 


0.0000 


C IUI 


m A 3/s 








0.0000 


D Ph(U) deg 




sameas 0 


Gas type 


0.0000 


E Hdot W 




ideal 


Solid type 


0.0000 


F Work W 




i 


27 — 










HXFRST 


Ambient HX 










sameas 3a 


a Area m A 2 


23.02 




A Ipl Pa 


0.6290 


b GasA/A 


103.2 




B Ph(p) 


deg 


5.0800E-03 


c Length m 


6.6930E-04 


C IUI m A 3/s 


1.7018E-03 


d yO m 


38.94 




DPh(U) 


deg 


-1.8309E-03 


e Heatln W G 


1.5886E-03 


EHdot 


W 


293.0 


fEst-T K (t) 


3.3455E-03 


F Work 


W 


sameas 0 


Gas type 


-1.8309E-03 


G Heat 


W 


copper 


Solid type 


293.0 




H MetalT 


K 


i 


28 — 










STKSLAB 


Prime Mover Stack 










2.6300E-03 


a Area m A 2 


16.26 




A Ipl Pa 


0.6220 


b GasA/A 


-43.67 




B Ph(p) 


deg 


2.4000E-02 


c Length m 


6.7080E-04 


C IUI m A 3/s 


3.0480E-04 


d yO m 


38.51 




D Ph(U) 


deg 


7.5000E-05 


e Lplate m 


1.5886E-03 


EHdot 


W 






7.4220E-04 


F Work 


W 
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293.0 


G T-beg K 


sameas 0 


Gas type 




293.0 


H T-end K 


kapton 

i 


Solid type 


.... 29 


-2.6033 


E-03 I StkWrk 


HXLAST 


Cold HX 








sameas 3a 


a Area 


m A 2 


16.63 


A Ipl Pa 


0.7050 


b GasA/A 




-43.90 


B Ph(p) deg 


3.0480E-04 


c Length 


m 


6.7062E-04 


C IUI m A 3/s 


1.2540E-03 


d yO 


m 


38.51 


D Ph(U) deg 


0.0000 


e Heatln 


W (t) 


7.3668E-04 


E Hdot W 


293.0 


fEst-T 


K =29H? 


7.3668E-04 


F Work W 


sameas 0 


Gas type 




-8.5I95E-04 


G Heat W 


copper 


Solid type 




293.0 


H MetalT K 



30 



SOFTEND 



0.0000 


a Re(Z) 


(t) 


0.0000 


A Ipl Pa 


0.0000 


b Im(Z) 


(t) 


0.0000 


B Ph(p) deg 








0.0000 


C IUI m A 3/s 








0.0000 


D Ph(U) deg 








0.0000 


E Hdot W 








0.0000 


F Work W 








0.0000 


G Re(Z) 


sameas 0 


Gas type 




0.0000 


H Im(Z) 


ideal 

t 


Solid type 


31 


0.0000 


I T K 


DIFFTAR 


Ipl mismatch 








0.0000 


a TargDi 


=31 A? 


8.0065E-05 A D1-D2 



1A b DIAddr 
30A c D2Addr 

! 32 

DIFFTAR Ph(p) mismatch 

0.0000 aTargDi =32A? 3.3780E-04 A D1-D2 

IB b DIAddr 
30B c D2Addr 

! 33 

DIFFTAR IUI mismatch 
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=33 A? 1.3884E-09 A D1-D2 



0.0000 a TargDi 

1C b DIAddr 
30C c D2Addr 

! 34 

DIFFTAR Ph(U) mismatch 

0.0000 a TargDi =34 A? -1.0691E-04 A D1-D2 

ID b DIAddr 
30D c D2Addr 

! The restart information below was generated by a previous run 
! You may wish to delete this information before starting a run 
! where you will (interactively) specify a different iteration 
! mode. Edit this table only if you really know your model! 
INVARS 504050607 27 5 
TARGS 5 29 6 31 1 32 1 33 1 34 1 
SPECIALS 0 
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APPENDIX B. THE MATHEMATICA PROGRAM FOR THE LEAST 
SQUARE FIT TO THE FREQUENCY RESPONSE 



(* This mathematica program read the plotting data file generated by *) 

(* DeltaE and perform a least square fit to find the resonance frequency *) 

Clear All["GlobaT*"] 

(* Read the data file generated by the DeltaE program. *) 

(* In this example, the data file is flowDT100.de *) 

datfile=ReadList[''flowDT100.de",Table[Number,{ 8 }]];(* There are eight columns in *) 

(* the data file *) 

(* Select the data range used in the least square fit, column 1 is the frequencies *) 

(* and column 7 is the pressure amplitude *) 

data=datfile[[Range[2,19],{ 1,7}]] (* Determine the bandwidth of fitting *) 

(* Initial guesses of A,f 0 , and Q *) 
soln = {A->150,f 0 ->437,Q->100}; 

Attributes[y]= { Listable } ; 

y[f_,A_,ffc_,Q_] := a (f/f 0 ) A 2 Abs[Cot[Pi (f/f 0 -I/(2 Q))]/(Pi (f/f 0 -I/(2 Q)))]; 

Attributes [yy] = {Listable}; 
yy[f_] = y[f,A,f 0 ,Q] /.soln; 
xlist = Transpose[data][[l]]; 
ylist = Transpose[data][[2]]; 

rmserror[A_,f 0 _,Q_] := Sqrt[ Apply[Plus,(ylist - y [xlist, A,f 0 ,Q]) A 2]]; 
rmserror[A,f 0 ,Q]/.soln; 

FindMinimum[rmserror[A,f 0 ,Q] , { A, { 1 00,200 } } , { f 0 , { 435,439 } } , { Q, { 90, 1 1 0 } } ] 

(* A is an arbitrary constant, f g the resonance frequency, Q the quality factor *) 

(* End of program *) 
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APPENDIX C. THE MATLAB PROGRAM FOR THE ANNULAR 
THERMOACOUSTIC PRIME MOVER 



Function updateF.m: The function performs the actual iteration to find the solutions for 
the constricted annular prime mover. 

% Begin of function % 

function [y,iterations]=updateF(y); 



% usage [newy, iterations] = updateF(y) 

% This function implements picards method for updating y 
% solving f(y,x)=0 for an annular resonator 

% The input y vector is a 8 by 1 vector: [Tm;pre;pim;Ure;Uim;fre;fim;H2dot] 

% only the last 5 elements (i.e. Ure, Uim, fre, fim and H2) are updated 

% newy = y - inv(J)*f(y,x); 

% where J is the Jacobian of f(y) 

% REQUIRES: fdjacF.m, funcvF.m, computeL.m, 

% Note: Before run this program, make sure to adjust the area ratio and select the 
% right temperature (Just active the line for the temperature profile) 

% Declare the global Variables 

global T2 T3 T4 T5 T6 T7 T8 T9 T10 T1 1 Th Tc Tmatrix Rs La Las R Rst r 
r = 2.565e-2; % Equivalent radius of the annular resonator, m 

Rs = 0.7; % The area ratio for the constriction 

% Compute the equivalent inductance of constriction 



d = 0.0529; 
h = 0.0497; 
dc = d*sqrt(Rs); 
he = h*sqrt(Rs); 
dst = 4.28e-2; 
hst = 6.12e-4; 

La = computeL(d,h,dc,hc); 
Las = computeL(d,h,dst,hst); 

R = sqrt((dc*hc)./(d.*h)); 



% Width of annulus, m 
% Height of annulus, m 
% Width of the constriction, m 
% Height of annulus, m 
% Width of a slit in the stack m 
% Height of a slit in the stack m 
% Compute the equivalent inductance of constriction 
% Compute the equivalent inductance of a slit in the 
% stack 
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Rst = sqrt((dst*hst)./(d.*h»; % Compute the equivalent area ratio for one slit 
% The measured temperature profiles for all the measurements 
% When run this main program, select the profile based on the porosity 
% The low mode and high mode share the same temperature profile 
% These temperatures are in °C 

% Low mode and high mode with constriction Rs=0.1 % Heater voltage 

% Tmatrix=[21 .0 21 21 21 21 21 21 21 21 20.6]; % heater = 0V 

% Tmatrix=[66.9 23.4 22.5 58.8 22.2 21.9 22.2 32.5 25.7 21.4]; % heater = 8V 

% Tmatrix=[ 100.7 25.8 24.3 86.8 23.8 23.2 23.7 41.2 30.4 22.6]; % heater = 11V 

% Tmatrix=[ 137.8 28.5 26.2 119.2 25.4 24.4 25.2 51.7 36.2 23.6]; % heater = 14V 

% Tmatrix=[218.6 36.6 32.4 194.5 30.6 28.7 30.2 82.1 51.8 27.3]; % heater = 20V 

% Tmatrix=[278.6 43.9 38 253.1 35.2 32.3 34.5 113.1 66.2 30.2]; % heater = 25V 
% Low Mode and high mode with constriction Rs=0.3 

% Tmatrix=[29.7 29.5 29.4 29.6 29.5 29.4 29.4 29.3 29.5 29.5]; % heater = 0V 

% Tmatrix=[74.3 30 29.1 66.8 28.9 28.5 28.8 37.8 32.2 28.4]; % heater = 8V 

% Tmatrix=[ 106.6 31.1 29.6 93.8 29.1 28.4 29.0 44.5 35.1 28.2]; % heater = 11V, 

% Tmatrix=[143.9 33.1 30.7 126 29.9 28.9 29.7 54.2 39.4 28.4]; % heater = 14V, 

% Tmatrix=[220.7 38.5 34.4 196.2 32.6 30.6 32.1 81.4 50.8 29.5]; % heater = 20V 

% Tmatrix=[278.4 44.9 39.1 251.5 36.3 33.2 35.5 110.3 63.1 31.5];% heater = 25V 

% Low and high Mode with constriction Rs=0.7 

% Tmatrix=[23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.5 23.4]; % heater = 0V 

% Tmatrix=[74.9 25.3 24.3 65.7 24 23.6 23.9 55.1 23.7 23.4]; % heater = 8V 

% Tmatrix=[ 108.5 27.2 25.5 94 24.9 24.2 24.8 74.5 30.4 23.9]; % heater = 11V 

% Tmatrix=[ 147.7 29.8 27.3 128.2 26.3 25.1 26.0 100.1 34.5 24.6];% heater = 14V 
% Tmatrix=[221.9 35.6 31.2 196 29.4 27 28.8 157.1 43.6 26.2]; % heater = 20V 

% Tmatrix=[280.7 42.5 36.2 252.9 33.3 29.7 32.4 219.2 53.5 29.5];% heater = 25V 
Tmatrix = Tmatrix+273; % Change temperature to Kelvin 
T2=Tmatrix(l); 

T3=Tmatrix(2); 

T4=Tmatrix(3); 

T5=Tmatrix(4); 

T6=Tmatrix(5); 

T7=Tmatrix(6); 

T8=Tmatrix(7); 

T9=Tmatrix(8); 

T10=Tmatrix(9); 
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T1 l=Tmatrix(10); 
Th = (2*T2+T5)/3; 
Tc = (2*T3+T4)/3; 
Nmax - 100; 

y=y0); 
y(l,l) = Th; 



% Weighted hot end targeted temperature, K 
% Weighted cold end temperature, K 
% Maximum number of integration 
% Make sure the input vector y is a column vector 
% Fixed the starting temperature at T h 



iterations=0; 

% Main loop for iterations 
while(l) 



iterations = iterations + 1 ; 






fvec=funcvF(y); 


% Compute the difference between the targets and % 




% calculated results 


% 


cheeky = norm(fvec./y(l:5)) 


% Check the norm of the ratio of the difference % 


if checky<1.0E-3 


% Set the converging threshold 


% 


fprintf(l,'Done\n'); 

y = y' 


% Print out the solution % 





[fvec,yend] = funcvF(y); 

yend = yend' % print out the integration results at the ends % 

Q = yend(l,6)/(2*yend(l,7)) % Find the quality factor % 
elseif iterations == Nmax % Terminate the iteration if iteration > N max % 
fprintf( 1 ,' 'Exceed max # of iterations'); 
return 



else 

J = fdjacF(y,fvec); % calculate the Jacobian % 
diffy = -JYfvec; % Use the LU factorization % 

if checky<1.0e-l 

diffy = diffy/2; % Close to converge, adjust the step size % 
elseif checky<5.0e-2 
diffy = diffy/4; 
elseif checky<1.0e-2 
diffy = diffy/6; 

end 

newy=y 

newy(4:8) = y(4:8)+diffy; % Perform a Newton iteration to update y % 

y=newy; % only update Ure, Uim, fre, fim and H2 % 

fprintf(l, 'Computing Update, iterations %4.0f \n',iterations) 
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end 



end 

% 

% End of updateF.m % 
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Function funcvF.m: The function integrate the coupled ordinary differential equations of 
the annular prime mover. 



% Begin of function % 

function [fvec,yend] = funcvF(yin) 
% 



% Includes the ends effect from high order mode method 

% This function integrates the annular resonator with the stack/heat exchanger assembly 
% at the end of the resonator (including ambient and hot heat exchanger) 

% constriction locates between xconl to xcon2 
% Requires: ductfunc2.m, ductfunc21.m, ductcon.m, 

% , stackfun.m, ode45.m, ambhxfunc.m, hothxfunc, zstuff.m 

% Declare the global variables 

global T2 T3 T4 T5 T6 T7 T8 T9 T10 T1 1 Th Tc Tmatrix Rs La Las R Rst 
yin=yin(:); 



d = 0.0529; % Width of annulus m 

h = 0.0497; % Height of annulus m 

dc = d*sqrt(Rs); % Width of the constriction m 

he = h*sqrt(Rs); % Height of constriction m 

cp=1005; % Isobaric heat capacity per unit mass 

% These are the positions of the thermocouples, m 



xp0= 0; 
xpl 1=0.00293; 
xpl2=0. 084345; 
xpl3=0. 39615; 
xpl4=0. 707955; 
xpl5=0.7556; 
xpl =0.76292; 
xp2 = 0.768; 
xp3 = 0.7920; 
xconl =0.1 33863; 
xcon2= 0.23290; 

xpL=0.7923; % Effective length of the annular resonator, m 

% These are not the resolution of the integration 
% These are the points in which the results are extracted 



% Thermocouple # 9 
% Thermocouple # 10 
% Thermocouple # 1 1 
% Thermocouple # 7 
% Thermocouple # 8 
% Ambient heat exchanger 
% Prime mover stack 
% Hot heat exchanger 

% 45° long constriction, 90° from the center of the stack assembly 
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xnl = 20; xn2 = 30; 

xn3 = 50; xn4 = 10; 

fvec = zeros(5,l); % Set up the fvec vector 

xllspan = linspace(xpO,xpl l,xnl); 

xl2span = linspace(xpl I,xpl2,xn2); 

xl3aspan = linspace(xpl2,xconl,xn2); 

xl3bspan = linspace(xconl,xcon2,xn2); 

xl3cspan = linspace(xcon2,xpl3,xn2); 

xl4span = Iinspace(xpl3,xpl4,xn3); 

xl5span = Iinspace(xpl4,xpl5,xn2); 

xl6span = linspace(xpl5,xpl,xnl); 

x2span = linspace(xpl,xp2,xnl); 

x3span = Iinspace(xp2,xp3,xn2); 

x4span = linspace(xp3,xpL,xnl); 

yO = yin(l:7); % make sure yO is column vector and only the first 7 elements used in a duct 
yO = y0(:); % Set up the starting boundary conditions 

[xll.yll] = ode45('ductfunc2 , ,xllspan,y0); % Use the MATLAB built in ODE45 function 
yOll = y 1 1 (xnl ,:); % Set up the boundary condition for next integration 

[x 1 2,y 1 2] = ode45('ductfunc21',xl2span,y01 1); 

y012 = yl2(xn2,:); % Set up the boundary condition for next integration 

[xl3a,yl3a] = ode45('ductfunc22al\xl3aspan,y012); 

% From duct to a constriction, U is continuous but p 2 =p,-U,Z a 

gamma =1.4; % Ratio, isobaric to isochoric specific heats 

rh = 1 .2825e-2; % Hydraulic radius of the resonator, rh=r/2, m 

rhc = rh*sqrt(Rs); % Hydraulic radius in the constriction, m 

omega = 2*pi*(yl3a(xn2,6) + j.*yl3a(xn2,7));% Extract the complex frequency 

rho = 353.065./yl3a(xn2,l); % Density of air as a function of temperature (K) 

c = 20.0447.*sqrt(yl3a(xn2,l)); % Speed of sound in air as a function of temperature (K) 

mue = 1 .846e-5 .*(y 1 3a(xn2, 1 )/300). A 1 .5 .*(4 1 0.4. /(y 1 3a(xn2, 1 )+l 1 0.4)); % Viscosity of air 

K=2.624e2.*(yl 3a(xn2,l)./300). A l .5. *(523.83 1 ./(yl3a(xn2,l)+... 

245.4. *exp(27.6./yl3a(xn2,l)))); % Thermal conductivity, Air 

Pr = 0.60928+0.2301 7. *exp(-0.0028565.*yl3a(xn2,l)); % Prandtl number 
dV = sqrt(2.*mue./(rho.*omega)); % Viscous penetration depth % 
dK = sqrt(2.*K./(rho.*cp.*omega)); % Thermal penetration depth % 
fv = (l-j)*dV/(2*rh); 
fvc = (l-j)*dV/(2*rhc); 
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fk = fv/sqrt(Pr); 
fkc = fvc/sqrt(Pr); 

k = omega/c *sqrt((l+(gamma-l)*fk)/(l-fv)); % Wave number in a regular duct 
kc = omega/c *sqrt((l+(gamma-l)*fkc)/(l-fvc)); % Wave number in the constriction 
[H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

Yl=diag(-j*sqrt(mymz2-k. A 2)./(k.*rho*c)); 

Y2=diag(-j*sqrt(nynz2-kc. A 2)./(kc.*rho*c)); 

Zimag=(HmO.'*((H*Y2*H.’+Yl)\HmO))/(Rs*4*r A 2); % Compute the acoustic impedance 

% caused by the constriction, using 
% higher order mode method 

Ral = rho.*omega.*dV.*(l-R)./R.*(l+(l-R. A 2)./(pi*R).*log((l+R)./(l-R)))/(Rs*4*r A 2); 
Zal = Ral+ Zimag; 

%Zal=0; % If do not want to included end correction, just set the impedance to zero 

% by activating this line 

U1 = yl3a(xn2,4) + j.*y 13a(xn2,5); % Volume velocity at entrance of the constriction 
y013b = yl3a(xn2,:); 

yOl 3b( 1 ,2) = yOl 3b( 1 ,2) - real(Zal*Ul); % Correct real part of pressure 

y013b(l,3) = yO 1 3b( 1 ,3) - imag(Zal*Ul); % Correct imag part of pressure 

[x 1 3b,y 1 3b] = ode45('ductcon',xl3bspan,y013b); % Constriction region 
% From a constriction to a duct , U is continuous but p 3 =p 2 -U 2 Z a 
rho = 353.065./y 1 3b(xn2, 1 ); 
c = 20.0447. *sqrt(yl3b(xn2,l)); 

mue = 1. 846e-5.*(yl3b(xn2,l)/300). A 1.5. *(410. 4./(yl3b(xn2,l)+l 10.4)); 

K = 2.624e-2.*(yl3b(xn2,l)./300) A 1.5.*(523.8306./(yl3b(xn2,l)+... 

245.4. *exp(-.27.6./yl3b(xn2,l))));% Thermal conductivity, Air 

Pr = 0.60928+0. 2301 7. *exp(-0.0028565.*yl3b(xn2,l));%Prandtl number 

dV = sqrt(2.*mue./(rho.*omega)); % Viscous penetration depth % 

dK = sqrt(2.*K./(rho.*cp.*omega)); % Thermal penetration depth % 

fv = (l-j)*dV/(2*rh); 

fvc = (l-j)*dV/(2*rhc); 

fk = fv/sqrt(Pr); 

fkc = fvc/sqrt(Pr); 

k = omega/c*sqrt(( l+(gamma- 1 )*fk)/( 1 -fv)); % Wave number in regular duct 
kc = omega/c*sqrt((l+(gamma-l)*fkc)/(l-fvc)); % Wave number in the constriction 
Yl=diag(-j*sqrt(mymz2-k. A 2)./(k.*rho*c)); 
Y2=diag(-j*sqrt(nynz2-kc. A 2)./(kc.*rho*c)); 
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Zimag=(HmO.'*((H*Y2*H.'+Yl)\HmO))/(Rs*4*r A 2); 

Ra2 = rho.*omega.*dV.*(l-R)./R.*(l+(l-R. A 2)./(pi*R).*log((l+R)./(l-R)))/(Rs*4*r A 2); 

Za2 = Ra2 + Zimag; 

%Za2 = 0; 

U2 = yl3b(xn2,4) + j*yl3b(xn2,5);% Volume velocity at start of the constriction 
y013c = yl3b(xn2,:); 

y013c(l,2) = y013c(l,2)-real(Za2*U2); % Correct real part of pressure 

y013c(l,3) = y013c(l,3)-imag(Za2*U2); % Correct imag part of pressure 

[xl3c,yl3c] = ode45('ductfunc22a2',xl3cspan,y013c); 
y013 = yl3c(xn2,:); 

[xl4,yl4] = ode45('ductfunc23',xl4span,y013); 
y014 = yl4(xn3,:); 

[xl5,yl5] = ode45(’ductfunc24',xl5span,y014); 
y015 = yl5(xn2,:); 

[xl6,yl6] = ode45('ductfunc25 l ,xl6span,y015); 
y016 = y 1 6(xn 1 

[x2,y2] = ode45(’ambhxfunc’,x2span,y016); % Integrate in the ambient heat exchanger 
% Compute the end corrections for the stack 

% Activate these command to include the end corrections in the stack region 
%nslit = 57; % Number of the slit 

%Aslit = 2.61936e-5; % Area for one slit, m 2 
%y0=3.06e-4; % Half spacing of the stack, m 

%ds = 4.28e-2; % Width of one slit in the stack region, m 

%hs = 6.12e-4; % Height of one slit in the stack region, m 

%omega = 2*pi*(y2(xnl,6) + j.*y2(xnl,7)); % Complex frequency 
%rho = 353.065 ./y2(xnl,l); 

%c = 20.0447. *sqrt(y2(xnl,l)>; 

%mue = 1.846e-5.*(y2(xnl,l)/300). A 1.5.*(410.4./(y2(xnl,l)+l 10.4)); 

%dV = sqrt(2.*mue./(rho.*omega)); % Viscous penetration depth in regular duct 
%dK = sqrt(2.*K./(rho.*cp.*omega)); % Thermal penetration depth in regular duct 
%fv = (l-j)*dV/(2*rh); 

%fk = fv/sqrt(Pr); 

%fvs = tanh((l+j)*yO/dV)/((l+j)*yO/dV); % fv for Parallel plate 

%fks = fvs/sqrt(Pr); % fk for Parallel plate 

%k = omega/c *sqrt((l+(gamma-l)*fk)/(l-fv)); % Wave number in regular duct 

%ks = omega/c *sqrt((l+(gamma-l)*fks)/(l-fvs));% Wave number in the constriction 
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%[H,HmO,mymz2,nynz2]=zstuff(d,h,ds,hs); 

%Yl=diag(-j*sqrt(mymz2-k. A 2)./(k.*rho*c)); 

%Y2=diag(-j*sqrt(nynz2-ks. A 2)./(ks.*rho*c)); 

%Zimag=(HmO.'*((H* Y2*H.'+Y 1 )\HmO))/(Aslit); 

%Rast3 = rho.*omega.*dV.*(l-Rst)./Rst.*(l+(l-Rst. A 2)./(pi*Rst).*log((l+Rst)./(l-... 
Rst)))/(Aslit); 

%Zastl = (Rast3 + Zimag)/nslit; % Effective impedance for end corrections of the 

% stack 

%Ustl = y2(xnl,4) + j.*y2(xnl,5); % Volume velocity at entrance of the constriction 
Zastl=0; % No end corrections for stack if this line is active 

y02 = [y2(xnl,:) yin(8)]; % Add the 8th element (i.e. H2) in the stack region 

%y02(l,2) = y02(l,2) - real(Zastl*Ustl); % Correct real part of pressure 

%y02(l,3) = y02(l,3) - imag(Zastl*Ustl); % Correct imag part of pressure 

[x3,y3] = ode45('stackfun',x3span,y02); % Integrate in the prime mover stack region 
% Activate these commands to include the end corrections in the stack region 
%Aslit = 2.61936e-5; % Area for one slit 

%rho = 353.065./y3(xn2,l); 

%c = 20.0447. *sqrt(y3(xn2,l)); 

%mue = 1 .846e-5.*(y3(xn2,l)/300). A l .5.*(410.4./(y3(xn2,l)+l 10.4)); 

%dV = sqrt(2.*mue./(rho.*omega)); % Viscous penetration depth in regular duct 
%dK = sqrt(2.*K./(rho.*cp.*omega));% Thermal penetration depth in regular duct 
%fv = (l-j)*dV/(2*rh); 

%fk = fv/sqrt(Pr); 

%fvs = tanh((l+j)*yO/dV)/((l+j)*yO/dV); % fv for Parallel plate 

%fks = fvs/sqrt(Pr); % fk for Parallel plate 

%k = omega/c*sqrt((l+(gamma-l)*fk)/(l-fv));% Wave number in regular duct 

%ks = omega/c*sqrt((l+(gamma-l)*fks)/(l-fvs));% Wave number in the constriction 

%[H,Hm0,mymz2,nynz2]=zstuff(d,h,ds,hs); 

%Yl=diag( -j*sqrt(mymz2-k. A 2)./(k.*rho*c)); 

%Y2=diag( -j*sqrt(nynz2-ks. A 2)./(ks.*rho*c)); 

%Zimag=(HmO.'*((H*Y2*H.'+Y 1 )\HmO))/(Aslit); 

%Rast4 = rho.*omega.*dV.*(l-Rst)./Rst.*(l+(l-Rst. A 2)./(pi*Rst).*log((l+Rst)./(l-... 
Rst)))/(Aslit); 

%Zast2 = (Rast4 + Zimag)/nslit;% Effective impedance for end corrections of stack 
Zast2=0; % If this line ic active, no end corrections for the stack is included 

%Ust2 = y3(xn2,4) + j.*y3(xn2,5); % Volume velocity at start of the constriction 
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y03 = v3(xn2,l:7): % only the first 7 elements are used in the hot heat exchanger 

%y03(1.2) = > 03(1.2) - reaI(Zast2*Ust2); 7c Correct real part of pressure 

%y03(l,3) = y03(l,3) - imag(Zast2*Ust2); 7c Correct imag part of pressure 

[x4,y4] = ode45(’hothxfunc\x4span,y03); 

vend = [y4(xnl,:> vin(8)]’; 7c This is the result at the end of the hot heat exchanger 
fvec(l.l) = vend(l,l)-Th; 7c Find the difference of the targeted hot temperature 
fvec(2:5,l)= yin(2:5,l)-yend(2:5,l);% Find the differences of the target pre, pim, Ure , Uim 
yend = vend(:): 7c Show results of the integration at the end 

% 

7c End of funcvF.m 
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Function fdjacF.m: This function computes the Jacobian using finite difference. 



% Begin of function % 

function df=fdjacF(yin,fin) 

% USAGE: df=fdjacF(yin,fin) 

% 

% Inputs: y, funcvF input vector to compute Jacobian about 
% f evaluation of funcv at y 

% Outputs: df = matrix with df(i,j) = dfuncv(i),dy(j) 

% This function compute an numerical approximation to the Jacobian matrix 
% for the function funcvF 
% 

% INPUTS: yin in the original guess vector 
% fin is the original differences of the target 

% Requires: funcvF.m 
% 



yin=yin(:); 

fin=fin(:); 

h=lE-10.*abs(yin)+eps; % make y perturbation be IE- 10 of y's current value 
% Perturb with respect to U r , U it f r , f, and H 2 . i.e. y(4:8) 

% DON'T perturb with respect to Tm, pr or pi 
df=zeros(5,5); % Initialize the Jacobian matrix (5 by 5) 
for i=4:8 

% perturb the only ith element of y % 

ypert=yin; 

ypert(i)=h(i)+yin(i); 

% find f with the perturbed y % 
pert = funcvF(ypert); 

% Forward difference formula % 
df(:,i-3) = (fpert(:)-fin(:))./h(i); % Jacobian array 
fprintf( 1 
end 
% 

% End of funcvF.m % 
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Function ductfunc2.m: This function contains the differential equations in a duct region. 



% Begin of function % 

function ydot = ductfunc2(x,y) 

% This function contains the differential equation in a dust region 
% dp,/dx = -jw*rho*Ul/((l-f v )*Area) 

% dU,/dx = -jw*Area*(l+(gamma-l)*f k )*pl/(rho*c A 2), U, is the volumetric velocity 
% The input y is a column vector with the component of 
% y [T m >Pre>Pim>!“^re>l“lj m TrcTim’F^] i 

% in a duct region, H 2 is not used, only y(l:7) is used in duct region 
% 

% 



global T2 T3 T4 T5 T6 T7 T8 T9 T10 T1 1 Th Tc r % All input temperature in C 
xpll = 0.00293; % Location of the thermocouple #9 

y(l) = Th; 
y=y(l:7); 

r_h = r/2 ; % Hydraulic radius of the resonator 

Area = 4*r A 2; % Cross-section area of the duct, m 2 



Tm = y(l); 
prl = y(2); 
pirn = y(3); 

Url = y(4); 

Uim = y(5); 
frl = y(6); 
fim = y(7); 

U = Url + j.*Uim; 
p = prl + j.*pim; 
f = frl + j*fim; 
w = 2*pi*f; 



% The volumetric velocity, s/m 3 
% Acoustic pressure 
% Complex frequency 



% Gas constants for air , Tm in Kelvin, at 1 ATM % 

% Thermal conductivity % 

K = 2.624e-2*(Tm/300) A 1.5*(523.8306/(Tm+245.4*exp(-27.6/Tm))); 

% Viscosity of Air % 

mue = 1.846e-5.*(Tm/300). A 1.5.*(410.4./(Tm+110.4)); 
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gamma =1.4; % Ratio, isobaric to isochoric specific heats 

cp=1005; % Isobaric heat capacity per unit mass 

% Compute gas properties and the f function 
[fv,fk,c,rho,K,Ksolid,Prandtl]=airproperty(Tm,w,l); 

Tmdot = (T9-Th)/xp 1 1 ; % Set the temperature gradient along the duct 

% set up the ODE equation in a duct region 

tempi = -j.*w.*rho.*U./((l-fv)*Area); 

temp2 = -j.*w.*Area.*(l+(gam-l).*fk).*p./(rho.*c. A 2)+... 

((fk-fv)*Tmdot*U)/((l-fv)*(l-Prandtl)*Tm); 

% Set up the ordinary differential equations 
prdot = real(templ); 
pidot = imag(templ); 

Urdot = real(temp2); 

Uidot = imag(temp2); 
frldot = 0; 
fimdot = 0; 

ydot =[Tmdot;prdot;pidot;Urdot;Uidot;frldot;fimdot]; 

% 

% End of ductfunc2.m % 
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Function ambhxfunc.m: This function contains the differential equations in the ambient 
heat exchanger. 

% Begin of function % 

function ydot = ambhxfunc(x,y) 

% This function contains the differential equations in the ambient heat exchanger 

% dp,/dx = -jw*rho*Ul/((l-f v )*Area) 

% dU/dx = -jw*Area*(l+(gamma-l)*f k )*pl/(rho*c A 2), U, is the volumetric velocity 

% The input y is a column vector with the component of 
% y-[T n) ,p re ,p im ,U re ,U, m) f re ,f im ,H 2 ] , 

% in the ambient heat exchanger region, H 2 is not used 
% Note that only y(l:7) is used here 
% 



global T2 T3 T4 T5 T6 T7 T8 T9 T10 T1 1 Th Tc % All input temperature in C 



y = y(l:7); 



y=y(:); 

Aambgas = 1.6547e-03; 
yOamb = 1.7018e-03; 
y(D=T6; 

Tm = y(l); 
prl = y(2); 
pirn = y(3); 

Url = y(4); 

Uim = y(5); 
frl = y(6); 
fim = y(7); 



% Make sure yO is a column vector 
% Gas area, m 2 

% Half spacing of plates for ambient heat exchanger, m 
% Ambient heat HX at T 6 Kelvin 



U = Url + j.*Uim; 
p = prl + j.*pim; 
f = frl + j*fim; 
w = 2*pi*f; 



% The volumetric velocity, s/m 3 
% Acoustic pressure 
% Complex frequency 



% Gas properties for air , Tm in Kelvin, at 1 atm % 

% Thermal conductivity 

K = 2.624e-2*(Tm/300) A l .5*(523.8306/(Tm+245.4*exp(-27.6/Tm))); 
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c = 20.0447*sqrt(Tm); % speed of sound in the air 

mue = 1.846e-5.*(Tm/300). A 1.5.*(410.4./(Tm+l 10.4)); % Viscosity of Air 

rho = 353.065/Tm; % Density of air 

gam = 1 .4; % Ratio, isobaric to isochoric specific heats 

cp=1005; % Isobaric heat capacity per unit mass 

deltaK = sqrt(2.*K./(rho.*cp.*w)); % Thermal penetration depth, m 
deltaV = sqrt(2.*mue./(rho.*w)); % Viscous penetration depth, m 

rootPrandtl=sqrt(mue.*cp./K); 

fv = tanh((l+j)*yOamb/deltaV)/((l+j)*yOamb/deltaV);% fv for Parallel plate 

fk = fv./rootPrandtl; % fk for Parallel plate 

% set up the ODE equations in the ambient heat exchanger region 

tempi = -j.*w.*rho.*U./((l-fv)*Aambgas); 

temp2 = -j.*w.*Aambgas.*(l+(gam-l).*fk),*p./(rho.*c. A 2); 

Tmdot = 0; % Note that it is assumed no temperature gradient across the ambient heat 

% exchanger 
prdot = real(templ); 
pidot = imag(templ); 

Urdot = real(temp2); 

Uidot = imag(temp2); 
frldot = 0; 
fimdot = 0; 

ydot =[Tmdot;prdot;pidot;Urdot;Uidot;frldot;fimdot]; 

% 

% End of ambhxfunc.m % 
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Function ductcon.m: This function contains the differential equations in the constricted 
duct. 

% Begin of function % 

function ydot = ductcon(x,y) 

% This function contains the differential equations in the constriction 
% dp,/dx = -jw*rho*Ul/((l-f v )*Area) 

% dU,/dx = -jw*Area*(l+(gamma-l)*f k )*pl/(rho*c A 2), U, is the volumetric velocity 
% The input y is a column vector with the component of 
^ y [T rn »Pre»Pim>hl re )Uj m ,f rc )fim»fl2] > 

% in a duct region, H2 is not used, only y(l:7) is used in duct region 

% 



global T2 T3 T4 T5 T6 T7 T8 T9 T10 T1 1 Th Tc Rs r 
% Define the locations of the constriction and thermoacoustics (in m) 
xpl2 = 0.084345; % Thermocouple # 10 

% Thermocouple # 1 1 
% Constriction locates at xconl to xcon2 
% 45° long constriction, 90° from center of stack 



xpl3 = 0.39615; 
xconl=0. 133863; 
xcon2= 0.23290; 



y=y( 1:7); 

r_h = r/2*sqrt(Rs) ; % Hydraulic radius of the constricted duct, m 

Area = Rs*4*r A 2; % Cross-section area of the duct, m 2 



Tm = y(l); 
prl = y(2); 
pirn = y(3); 

Url = y(4); 

Uim - y(5); 
frl = y(6); 
fim = y(7); 

U = Url + j.*Uim; 
p = prl + j.*pim; 
f = frl + j*fim; 
w = 2*pi*f; 



% Volumetric velocity, s/nv 
% Acoustic pressure 
% Complex frequency 



% Gas properties for air , Tm in Kelvin, at 1 atm % 



152 



cp=1005; % Isobaric heat capacity per unit mass 

c=20.0447.*sqrt(Tm); % Speed of sound in air 
rho=353.065./Tm; % Density of air 

gam =1.4; % Ratio, isobaric to isochoric specific heats 

% Thermal conductivity. Air, Tm in Kelvin 

K = 2.624e-2.*(Tm./300). A 1.5.*(523.8306./(Tm+245.4.*exp(-27.6./Tm))); 
mue = 1 ,846e-5.*(Tm/300). A l .5.*(410.4./(Tm+l 10.4)); % Viscosity of Air 

Prandtl = 0.60928+0.2301 7. *exp(-0.0028565.*Tm); % Prandtl number 



deltaK = sqrt(2.*K./(rho.*cp.*w)); 
deltaV = sqrt(2.*mue./(rho.*w)); 
fv = (l-j).*deltaV./(2*r_h); 
fk = fv./sqrt(Prandtl); 

Tmdot = (T 1 1 -T 1 0)/( xpl3-xpl2); 



% Thermal penetration depth in the constriction, m 
% Viscous penetration depth in the constriction, m 
% f v of the constriction 
% f k of the constriction 

% Temperature gradient across the constriction 



% set up the ODE equation in a constricted duct region 

tempi = -j.*w.*rho.*U./((l-fv)*Area); 

temp2 = -j.*w.*Area.*(l+(gam-l).*fk).*p./(rho.*c. A 2)+... 

((fk-fv)*Tmdot*U)/((l-fv)*(l-Prandtl)*Tm); 
prdot = real(templ); 
pidot = imag(templ); 

Urdot = real(temp2); 

Uidot = imag(temp2); 
frldot = 0; 
fimdot = 0; 



ydot =[Tmdot;prdot;pidot;Urdot;Uidot;frldot;fimdot]; 

% End of ductcon.m % 
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Function stackfun.m: This function contains the differential equations in the prime 

mover stack. 

% Begin of function % 

function ydot = stackfun(x,y) 

% This function contains the differential equations in the prime mover stack 
% to integrate in the stack region 

% The input y is a column vector with the component of 
% y=[T m ;p f e;p, m ;U ra ;U im ;f re ;f jn) ;H 2 ]; 

% 

% Requires: airproperty.m 
% 



global T2 T3 T4 T5 T6 T7 T8 T9 T10 T1 1 Th Tc 



y=y(0; 

Agas = 1.7995e-3; 
AsAgas = 0.2237; 
y(l) = Tc; 

Tm = y(l); 
prl = y(2); 
pim = y(3); 

Url = y(4); 

Uim = y(5); 
frl = y(6); 
fim = y(7); 

H2 = y(8); 

U = Url + j.*Uim; 
p = prl + j.*pim; 
f = frl + j*fim; 



% Total gas cross section area, m 2 , 56 cover glasses 
% Asolid/Agas in the stack region 
% Ambient side's temperature of the prime mover stack 



% Time-averaged energy flux, constant along the stack 
% Volumetric velocity, s/m 3 
% Acoustic pressure 
% Complex frequency 



w = 2*pi*f; 

% Gas properties for air 

[fv,fk,c,rho,K,Ksolid,Prandtl] = airproperty(Tm,w,2); % index == 2 for parallel plates 

% geometry 

gam =1.4; % Ratio, isobaric to isochoric specific heats 

cp=1005; % Isobaric heat capacity per unit mass 
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% H2 is the time-averaged energy flux, constant along the stack % 
dTmdxl = H2/Agas-0.5/Agas*real(p*U'*(l-((fk-fv’)/((l+Prandtl)*(l-fv'))))); 
dTmdx2 = 0.5 *rho*cp*abs(U) A 2*imag(fk+Prandtl*fv')/(real(w)*( 1 -Prandtl A 2) *abs( 1 . 
fv) A 2*Agas A 2) - K- AsAgas*Ksolid; 

dTmdx = dTmdxl/dTmdx2; % Temperature gradient in the stack region 
tempi = -j.*w.*rho.*U./((l-fv)*Agas); 

temp2 = -j.*w.*Agas.*(l+(gam-l).*fk).*p./(rho.*c. A 2)+((fk-fv)*dTmdx*U)/((l-fv)*(l 
Prandtl)*Tm); 

Tmxdot - dTmdx; 
prdot = real(templ); 
pidot = imag(templ); 

Urdot = real(temp2); 

Uidot = imag(temp2); 
frldot = 0; 
fimdot - 0; 

H2dot = 0; % H2 is constant through the stack 

ydot =[Tmxdot;prdot;pidot;Urdot;Uidot;frldot;fimdot;H2dot]; 

% End of stackfun.m % 
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Function airproperty.m: This function compute the properties for air. 

% Begin of function % 

function [fv,fk,c,rho,K,Ksolid, Prandtl] = airproperty(T,w,index) 

% This function compute the air properties at 1 atm and T Kelvin 

% if index == 1 » BL approximation for regular duct 

% index == 2 » Parallel plate for stack region 

% index == 3 » BL approximation for constricted duct 

global Rs 

r_h = 2.5654E-2/2; % Hydraulic radius of the resonator is r/2 

yO = 3.06e-4; % Half spacing of plates for stack, m 

cp=1005; % Isobaric heat capacity per unit mass 

c=20.0447.*sqrt(T); % Speed of sound in air 

rho=353.065./T; 

% Thermal conductivity, Air 

K = 2.624e-2.*(T./300). A 1.5.*(523.8306./(T+245.4.*exp(-27.6./T))); 

Ksolid = 0.2*(l-exp(-T/100)); % Thermal conductivity, kapton 
mue = 1 . 846e-5 . * (T/300) . A 1 .5 . *(4 1 0.4./(T + 1 1 0.4)) ; % Viscosity of Air 
Prandtl = 0.60928+0. 23017. *exp(-0. 0028565. *T); % Prandtl number 
deltaK = sqrt(2.*K./(rho.*cp.*w)); % Thermal penetration depth, m 

deltaV = sqrt(2.*mue./(rho.*w)); % Viscous penetration depth, m 

if index == 1 

fv = (l-j).*deltaV./(2*r_h); % fv for boundary layer approximation 

elseif index == 2 

fv = tanh((l+j)*yO/deltaV)/((l+j)*yO/deltaV); % fv for Parallel plate 
elseif index == 3 

fv = (l-j).*deltaV./(2*r_h*sqrt(Rs));% fv for boundary layer approximation 

% in the constricted duct 

else 

fprintf(l,' Wrong index number !!') 
return 

end 

fk = fv./sqrt(Prandtl); 

% End of airproperty.m % 
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Function plotpUF.m: This function plot the acoustic pressure distribution in the annular 
prime mover prime mover stack 

% Begin of function % 

function [x,y,pnorm]= plotpUF(yin) 

% This function plot the acoustic pressure amplitude of the auunlar prime mover 
% It uses the solution from “updataF.m” and integrates "yin" from x=0 to x=xL 
% 

% Note: The first section of this function which is not included here, 

% contains the function ’’funcvF.m”. It uses funcvF.m to integrate 

% the solution from “updateF.m” from xO to xL 



% 

% Put :funcvF.m here 
% 



dT = Th-Tc; % Temperature difference across the prime mover stack, K 

freq = y0(6,l) +j*y0(7,l); % Complex frequency 

% Combine the solution vector 
% x is the position vector, y is the solution matrix 
x=[xl I;xl2;xl3a;xl3b;xl3c;xl4;xl5;xl6;x2;x3;x4]; 

y=[yll(:,l:7);yl2(:,l:7);yl3a(:,l:7);yl3b(:,l:7);yl3c(:,l:7);yl4(:,l:7);yl5(:,l:7);... 

y 16(:,l:7);y2(:,l:7);y3(:,l:7);y4(:,l:7)]; 

T = y(:,l); 

p =y(:,2)+j*y(:,3); % Complex acoustic pressure 

U =y(:,4)+j*y(:,5); % Complex volume velocity 

pm=zeros(8,l); 

% Find the theoretical pressure amplitude at the locations of each microphones 
% There are 8 microphones used in the measurement 
for 1=1:8 

xm(i) = 0.03483+(i-l )*9.90375e-2; % Microphone location, m 

[Y, index] = min(abs(x-xm(i))); % Compute the corresponding index 

pc(i)=pnorm(index); % Compute the theoretical values at 

% each microphones’ locations 



end 
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pc=pc(:); 

% These are the calibrated measured pressure amplitude at each microphones 
% For the constricted annular prime mover 

% Note that all these measurements are subjected to the measured temperature profile listed 
% in the function “updateF.m” 



% 

% Low mode, With constriction 
% Rs = 0.7, low mode 

%pm=[34.64;20.63;14.82;31.92;31.18;13.42;12.07;28.41]'; 
%pm=[35.44;21.3;15;32.63;31.94;13.78;12.16;29.07]'; 
%pm=[35 .72 ;2 1 .77 ; 1 4.95 ;3 3.05 ;32.5 3; 1 4.09; 1 2.17 ;29.39]*; 
%pm=[36.26;22.37;15.14;33.77;33.16;14.5;12.3;30.03]'; 
%pm=[37.38;23.76;15.3;35.39;35.24;15.57;12.62;31.45]'; 
%pm=[37.98;24.63;15.46;36.71 ;36.71;16.29;13.1 ;32.65]'; 



% Heater = 0V 
% Heater = 8V 
% Heater = 1 1 V 
% Heater =14V 
% Heater = 20V 
% Heater = 25 



% Rs=0.3, low mode 

%pm=[34.89;26.09;12.32;34.67;28.18;10.54;11.6;26.91]'; 
%pm=[39.74;29.49;14.06;38.56;31.29;l 1.57;13.08;29.95]'; 
%pm=[43.1 1;32.27;15.82;42.16;34.18;12.56;14.41;32.77]'; 
%pm=[47.39;35.68;18.16;46.65;37.88;13.71 ;16.1 ;36.31]'; 
%pm=[58.26;45.57;24.5;60.3;48.71 ;17.44;21.35;47.11]'; 
%pm=[55.92;55.06;30.23;74.79;60.53;21.48;27.1;58.8]'; 



% Heater = 0V 
% Heater =V8 
% Heater = 11V 
% Heater = 14V 
% Heater = 20V 
% Heater = 25 V 



% Rs=0.1, low mode 

%pm=[40.57;37.74;34.13;38.99;28.57;10.43;9.76;26.8]'; 

%pm=[47.45;44.12;39.87;45.48;33.24;11.91;11.76;31.61]’; 

%pm=[53.45;50.17;45.87;52.2;38;13.43;13.95;36.4]'; 

%pm=[62.9;59.89;54.71;62.74;45.34;15.77;17.47;43.7]'; 

%pm=[91.08;97.23;89.03;104.3;75.09;25.35;32.76;73.7]'; 

%pm=[l 17.4; 155.5 ; 145.9 ; 178.6; 128.7 ;43.3;60.38; 135.1]'; 



% Heater = 0V 
% Heater = 8V 
% Heater = 11V 
% Heater = 14V 
% Heater = 20V 
% Heater = 25 V 



% High mode, With constriction 
% Rs=0.7, high mode 

%pm=[7.29;12.45;l 1.79;6.47;4.78;9.75;12.213;6.14]'; 
pm=[8.03;13.09;l 1 .48 ;5.6;5.44; 1 0. 1 9; 1 1 .67 ;5. 1 7]'; 
%pm=[8.14;12.99;l 1.19;5.43;5.66;10.14;1 1.3;4.92]'; 
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% Heater = 0V 
% Heater = 8V 
% Heater = 11V 



%pm=[8.41;13.22;l 1.19;5.39;5.97;10.3;11.23;4.84]'; 
%pm=[8.7;13.29;10.96;5.38;6.43;10.35;10.94;4.75]'; 
%pm=[8. 59; 1 3.07; 10. 67;5. 5;6.65; 10.26; 10.7;4.79]’; 
% Rs=0.3, high mode 

%pm=[6.86;12.83;6.18;6.49;2.87;9.72;11.99;4.47]’; 
%pm=[7.24;13.41;6.43;6.56;3.05;10;12.12;4.32]'; 
%pm=[7.38;13.64;6.69;6.62;3.19;10.18;12.21;4.22]’; 
%pm=[7.45;13.83;6.95;6.63;3.32;10.34;12.27;4.06]'; 
%pm=[7.41;14.4;7.48;6.81;3.67; 10.87; 12.7 1;3. 87]'; 
%pm=[5.56;14.5;7.61 ;6.89;3.89;11.16;12.93;3.67]'; 
% Rs=0.1, high mode 

%pm=[8.75;15.87;14.26;8.5;2.63;l 1.37;12.68;4.99]'; 

%pm=[9;16.27;14.48;8.59;2.83;11.67;12.85;4.76]’; 

%pm=[8.98;16.33;14.65;8.6;2.94;l 1.79;13;4.56]’; 

%pm=[8.98;16.8;15.07;8.84;3.13;12.21;13.66;4.48]'; 

%pm=[7.69;17.31;15.12;8.99;3.5;12.87;15.10;4.17]'; 

%pm=[6.39;17.4;14.87;9.46;3.77;13.43;15.79;4.01]'; 



% Heater = 14V 
% Heater = 20V 
% Heater = 25V 

% Heater = OV 
% Heater = 8V 
% Heater = 1 IV 
% Heater = 14V 
% Heater = 20V 
% Heater = 25V 

% Heater = OV 
% Heater = 8V 
% Heater = 11V 
% Heater = 14V 
% Heater = 20V 
% Heater = 25V 



% A least fit is applied to the predicted values and the measured data 
A=pm*pc/(pc'*pc); % Find the least square fit constant 
pmscal=pm./A; % Normalized microphone voltage by the constant A 

% These are the actual positions of the microphones, m 

xpm=[0. 03483;0.13387;0. 2329 1;0.331 95 ;0.43099;0. 53003 ;0. 62907 ;0. 7281 1]; 
clg 

lx =[0.7629 0.7629]; % Mark the starting location for stack assembly 

ly = [0.1.2]; 

lxl=[0. 133863 0.133863 ]; % Mark the location of the constriction 

ly 1 =[0 1.2]; 

lx2=[0.2329 0.2329 ]; % Mark the location of the constriction 

ly2=[0 1.2]; 

% Plot the predicted pressure amplitude vs. the least square of the measured data 
plot(x,abs(p)./max(abs(p)), , b- , ,xpm,pmscal, , ro',lx,ly,'g--’,lxl,lyl,'g:’,lx2,ly2,’g:') 
legendf'Matlab program ','Measured data',' Stack region', 'Constriction'); 

%title([’Low mode with end effect 3 ; dT = ’,num2str(dT),' Kelvin ;... % Title for low mode 
% freq = ’,num2str(freq),' ; Rs = ',num2str(Rs)]) 

title([’High mode with end effect 3 ; dT = ’,num2str(dT),' Kelvin ;... % Title for high mode 
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freq = ',num2str(freq),' ; Rs = 
xlabel('Position, m') 
ylabel('Normalized pressure') 
axis([0 0.7923 0 1.2]) 


- ',num2str(Rs)]) 


% End of plotpUF.m % 
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Function computeL.m: This function computes the equivalent inductance for the series 
acoustic impedance Za(w) of a size discontinuity between a duct dxh and dcxhc using 
higher order mode theory. 

This function is contributed by Ralph T. Muehleisen 

% Begin of function % 

function La=computeL(d,h,dc,hc); 

% 

% Inputs 

% d = depth of main channel 
% h = height of main channel 
% dc = depth of constriction 
% he = height of constriction 
% Outputs 

% La = equivalent inductance of constriction 
% 

% This function computes the equivalent inductance for the 
% series acoustic impedance Za(w) of a size discontinuity between a 
% duct dxh and dcxhc using higher order mode theory. 

% The discontinuity is assumed to be symmetric 
% in the y direction and asymmetric in the z direction 
% (i.e. going from alxbl to a2xb2 the open area in duct 2 is 
% (d-dc)/2<y<(d+dc)/2 and (h-hc)<z<h 
% 

% To use this add the following lines to your code 
% 

% In the beginning: 

% La = computeL(d,h,dc,hc) 

% R=sqrt((dc*hc)./(d.*h)); 

% Ra = rho.*w.*dV.*(l-R)./R.*( 1+ (l-R. A 2)./(pi*R).*Iog((l+R)./(l-R)) ); 

% Za = Ra + j*w*La; 

R=sqrt((dc*hc)./(d.*h)); 

% compute the coupling matrix H, using n2=3 modes in small constriction 



161 



% and select n 1 based on relative sizes 
n2=3; 

n 1 =round(n2/sqrt(R)); 



% call h2d to get the full coupling matrix 
[Hfull,rny,mz,ny,nz] = h2d(R,nl,n2); 

[MH,NH] = size(Hfull); 

H00= Hfull (1,1); 

% extract out the higher order mode terms 
H = Hfull(2:MH,2:NH); 

% Extract out the plane wave terms 

HmO = Hfull(2:MH,l); 

my = my(2:MH); 

mz = mz(2:MH); 

ny = ny(2:NH); 

nz = nz(2:NH); 

mymz2 = (my.*pi./d) A 2 + (mz.*pi./h). A 2; 
nynz2 = (ny.*pi./dc). A 2 + (nz.*pi./hc). A 2; 

% The inductance approximation, it assumes k2 « my2mz2 or ny2nz2 and k~kc 

Yin = diag(sqrt(mymz2)); 

Y2n = diag(sqrt(nynz2)); 

La = rho* (HmO.' *inv(H*Y2n*H.' + Yln)*HmO)/(hc*dc); 
return 

% End of computeL.m % 
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Function zstuff.m: This function computes the scattering matrices used in computation of 
series acoustic impedance Za(w) of a size discontinuity between a duct dxh and dcxhc 
using higher order mode theory. 

This function is contributed by Ralph T. Muehleisen 

% Begin of function % 

function [H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

% usage [H,HrnO,rnymz2,nynz2]=zstuff(d,h,dc,hc); 

% 

% Inputs 

% d = depth of regular duct 
% h = height of regular duct 
% dc = depth of constriction 
% he = height of constriction 
% 

% Outputs 

% H = higher order mode scattering matrix 
% HmO = scattering vector back to plane wave modes 
% mymz2, nynz2 = indices of higher order modes 
% 

% This function computes the scattering matrices used in computation 
% of series acoustic impedance Za(w) of a size discontinuity between a 
% duct dxh and dcxhc using higher order mode theory. 

% The discontinuity is assumed to be symmetric 
% in the y direction and asymmetric in the z direction 
% (i.e. going from alxbl to a2xb2 the open area in duct 2 is 

% d-dc)/2<y<(d+dc)/2 and (h-hc)<z<h 

% 

% Add the following lines to your program 
% somewhere in the beginning add 
% [H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

% 

% In the frequency dependent part add following. If the constriction is on the right 
% Yl=diag( -j*sqrt( mymz2-k. A 2)./(k*rho*c)); 

% Y2=diag( -j*sqrt(nynz2-kc. A 2)./(kc.*rho*c)); 
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% Zimag=(HmO.'*((H*Y2*H.’+Yl)\HmO))/Sc; 

% R=sqrt((dc*hc)./(d.*h)) 

% Ra = rho.*w.*dV.*(l-R)./R.*(l+(l -R. A 2)./(pi*R).*log(( 1+R)./(1-R))); 

% Za = Ra+ Zimag; 

% 

R=sqrt((dc*hc)./(d.*h)); 

% Compute the coupling matrix H, using n2=3 modes in small constriction 

% and select nl based on relative sizes 

n2=3; 

nl=round(n2/R); 
if nl> 10 

nl=10; 

end 

% Call h2d.m to get the full coupling matrix 
[Hfull,my,mz,ny,nz]=h2d(R,nl,n2); 

[MH,NH]=size(Hfull); 

H00=Hfull(l,l); 

% Extract out the higher order mode terms 
H=Hfull(2:MH,2:NH); 

% extract out the plane wave terms 

HmO =Hfull(2:MH,l); 

my=my(2:MH); 

mz=mz(2:MH); 

ny=ny(2:NH); 

nz=nz(2:NH); 

mymz2 = (my.*pi./d). A 2 + (mz.*pi./h). A 2; 
nynz2= (ny.*pi./dc). A 2 + (nz.*pi./hc). A 2; 
return 

% End of zstuff.m % 
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Function h2d.m: This function is used to calculate the coupling matrix H for a 
change of size of a square duct (axa)->(bxb) with R=b/a. 

This function is contributed by Ralph T. Muehleisen 

% Begin of function % 

function [H,my,mz,ny,nz]=h2d(R,N ] ,N2) 

% 

% Inputs: 

% Ry = a2/al = width ratio of duct, 

% Rz = b2/bl = height ratio of duct 

% Nl= # of modes in large duct region and 
% N2 = # of modes in small duct in one direction. 

% 

% Outputs: 

% H = coupling matrix 
% Mx = large duct y mode matrix 
% Mz= large duct z mode matrix 
% Ny = small duct y mode matrix 

% Ny = small duct z mode matrix 

% 

% The routine uses the N1 and N2 lowest frequency modes in each section. 

% It assumes that the y direction constriction is symmetric (i.e b is 
% centered on a) and the z constriction is asymmetric (b is at one 
% edge of a 
% 

Nt2=N2. A 2; 

Ntl=Nl. A 2; 

% Generate ny, nz 

temp=(0:N2-l)'; 

ny=zeros(Nt2,l); 

nz=ny; 

for i=0:N2-l 

ny(i*N2 + (temp+l),l)=[i.*ones(size(temp))]; 
nz(i*N2 + (temp+l),l)=temp; 
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end 



% now sort the modes by frequency 



f2=sqrt(ny. A 2 + nz. A 2); 

[y,I]=sort(f2); 

ny=ny(I); 



% Compute relative freq of each mode 
% Sort the freq. 

% Pull out ny and nz 



nz=nz(I); 

% Generate my,mz 

temp=(0:Nl-l)'; 

my=zeros(Ntl,l); 

mz=my; 

for i=0:Nl-l 

my(i*Nl + (temp+l),l)=[i.*ones(size(temp))]; 
mz(i*Nl + (temp+l),l)=temp; 



mz=mz(I); 

% convert my,mz,ny,nz into the matrices used to compute H=Hy*Hy; 
[Ny,My]=meshgrid(ny,my); 

[Nz,Mz]=meshgrid(nz,mz); 

Mypio2=My . *pi/2; 

Hy=zeros(size(My)); 

Hz=Hy; 

% first compute H for the symmetric y direction 

e=(l-R)/2; 

sme=sin(My*pi*e); 

cme=cos(My*pi*e); 

smepr=sin(My*pi*(e+R)); 

smep2r=sin(My*pi*(e+2*R)); 

% first fill in the general terms 

Hy=-2*R. A (3/2).*My.*(sme-(-l ). A Ny.*smepr)./(pi*((My.*R). A 2-Ny. A 2+eps»; 
% find the n=m*R terms (including m=n=0) and set them to (-l) A (m+n)*Sqrt(R) 
I=find(abs(Ny-R*My)<10*eps); 

temp=(2*My.*pi.*R.*cme-sme+smep2r)./(4*My*pi*sqrt(R)+eps); 



end 

fl=sqrt(my. A 2 + mz. A 2); 

[y,I]=sort(fl); 

my=my(I); 



% Compute relative freq of each mode 
% Sort the freq. 

% Pull out ny and nz 
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Hy(I)=temp(I); 

% find the m=0 terms and set them to zero 
I=find(My==0); 

Hy(I)=0*I; 

% Find the n=0 terms, they were computed above by the general expression 

% but are sqrt(2) too big 

I=find(Ny==0); 

Hy(I)=Hy(I)/sqrt(2); 

% find the m=n=0 terms and set them equal to sqrt(R) 

I=find(My==0 & Ny==0); 

Hy(I)=sqrt(R)*ones(size(I)); 

% second compute H for the asymmetric z direction 
e=(l-R); 

sme=sin(Mz*pi*e); 

cme=cos(Mz*pi*e); 

smepr=sin(Mz*pi*(e+R)); 

smep2r=sin(Mz*pi*(e+2*R)); 

% first fill in the general terms 

Hz=-2*R. A (3/2).*Mz.*(sme-(-l). A Nz.*smepr)./(pi*((Mz.*R). A 2-Nz. A 2+eps)); 

% find the n=m*R terms (including m=n=0) and set them to (-l) A (m+n)*Sqrt(R) 
I=find(abs(Nz-R*Mz)<10*eps); 

temp=(2*Mz.*pi.*R.*cme-sme+smep2r)./(4*Mz*pi*sqrt(R)+eps); 

Hz(I)=temp(I); 

% Find the m=0 terms and set them to zero 
I=find(Mz==0); 

Hz(I)=0*I; 

% find the n=0 terms, they were computed above by the general expression 

% but are sqrt(2) too big 

I=find(Nz==0); 

Hz(I)=Hz(I)/sqrt(2); 

% find the m=n=0 terms and set them equal to sqrt(R) 

I=find(Mz==0 & Nz==0); 

Hz(I) = sqrt(R)*ones(size(I)); 

H = Hy.*Hz; 

% End of h2d.m % 
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Function update2st.m: This function iterates the ordinary differential equation for a two- 
stack annular prime mover. 

% Begin of function % 

function [y,iterations]=update2st(y); 

% 

% 

% usage [newy,iterations] = update2st(y) 

% 

% 

% This function implements picards method for updating y 
% solving f(y,x)=0 for a resonator 

% The input y vector is a 8 by 1 vector. [T m ;p re ;p im ;U rc ;U ira ;f re ;f jm ;H 2a ;H 2b ] 

% only the last 6 elements (i.e. U re , U im , f re , f im , H 2a , and H 2b ) are updated 
% 

% newy = y - inv(J)*f(y,x); 

% where J is the Jacobian of f(y) 

% REQUIRES: fdjac2st.m, func2st.m 

Nmax = 50; 

y = y(0; 

iterations = 0; 
while(l) 

iterations = iterations + 1; 
fvec=func2st(y); 

cheeky = norm(fvec./y(l:6)) % check the norm of the ratio of the diffenence % 
if cheeky < 1.0E-2 % to the solution, be sure to use "./" not 7" % 

fprintf( 1 ,'Done\n'); 

y = y' % Print out the solution % 

[fvec,yend] = func2st(y); 

yend = yend' % print the integration results at the ends % 

Q = yend(l,6)/(2*yend(l,7)) % Find the quality factor % 
return 

elseif iterations == Nmax 

fprintf(l, 'Exceed max # of iterations'); 
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return 



else 

J = fdjac2st(y,fvec); % calculate the Jacobian % 

diffy = -J\fvec; % Use LU factorization to find the correction to the guess % 

if checkycl .Oe-1 % Reduce the step size when close to solution % 

diffy = diffy/2; 
elseif cheeky < 5.0e-2 
diffy = diffy/4; 
elseif cheeky < 5.0e-3 
diffy = diffy/6; 

end 

newy = y; 

newy(4:9) = y(4:9) + diffy; % Do a Newton iteration to update y % 
y = newy % only update Ure, Uim, fre, fim, H2a, and H2b % 

fprintf( 1 'Computing Update, iterations %4 .Of \n', iterations) 
end 

end 

% End of update2st.m % 
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Function fdjac2st.m: This function computes the Jacobian for the two-stack annular 
prime mover using finite difference. 

% Begin of function % 

function df=fdjac2st(yin,fin) 

% USAGE: df=fdjac2st(yin,fin) 

% 

% INPUTS: y funcv input vector to compute Jacobian about 
% f evaluation of funcv at y 
% OUTPUTS: df = matrix with df(i,j) = dfuncv(i),dy(j) 

% This function compute an numerical approximation to the Jacobian matrix 
% for the function funcv 
% INPUTS: 

% yin in the original guess vector 
% fin is the original differences of the target 

% 

% Requires: func2st.m 



yin = yin(:); 
fin = fin(:); 

h = lE-10.*abs(yin)+eps; % make y perturbation be IE-10 of y's current value 
% Perturb with respect to Ur, Ui, fr, fi and H2a H2b . i.e. y(4:9) 

% DON'T perturb with respect to Tm, pr or pi 

df = zeros(6,6); % Initialize the Jacobian matrix (6 by 6) 

for I = 4:9 

% perturb the only ith element of y 
ypert=yin; 

ypert(i) = h(i)+yin(i); 

% find f with the perturbed y 
fpert = func2st(ypert); 

% Forward difference formula % 

df(:,i-3) = (fpert(:)-fin(:))./h(i); % Jacobian matrix % 
fprintf(l,V); 

end 

% End of fdjac2st.m % 
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Function funcv2st.m: The function integrate the coupled ordinary differential equations 
of the two-stack annular prime mover. 



% Begin of function % 

function [fvec,yend] = func2st(yin) 

% This function integrates coupled ODE for the two stack annular prime move 
% Requires: 

% duct2stl.m, hothx2stl.m, stackl.m, stack2.m, ambhx2stl.m, ambhx2st2.m, 

% hothx2stl.m, hothx2st2.m, ode45.m 

% Note: Most of these functions are similar to the functions for the constricted 

% annular prime mover. Therefore, only the main functions are included here. 



global Th Tc % All input temperature in K 

% Before to execute the program, be sure to adjust % 

% the target hot heat exchanger temperature to the desired value % 
Th = 393; 



Tc = 293; 
yin = yin(:); 
xpO = 0; 
xpl = 0.16869; 
xp2 = 0.1689948; 
xp3 = 0.1929948; 
xp4 = 0.198075; 
xp5 = 0.76292; 
xp6 = 0.768; 
xp7 = 0.792; 
xpL = 0.7923; 
xnl = 50; 



% Starting location for the first satck assembly 

% End of the first stack assembly 

% Starting location for the second satck assembly 

% End of the second stack assembly 
% Effective length of the resonator, m 



xn2 = 10; 

xn3 = 20; 

fvec = zeros(6,l); 

xlspan = linspace(xpO,xpl,xnl); 

x2span = Iinspace(xpl,xp2,xn2); 

x3span = Iinspace(xp2,xp3,xn3); 

x4span = Iinspace(xp3,xp4,xn2); 
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x5span = Iinspace(xp4,xp5,xnl); 
x6span = Iinspace(xp5,xp6,xn2); 
x7span = Iinspace(xp6,xp7,xn3); 
x8span = Hnspace(xp7,xpL,xn2); 

% 

% Make sure that yO is column vector, only the first 7 elements are used in the duct% 
yO = yin( 1:7); 
yO = y0(:); 

[x 1 ,y 1] = ode45('duct2str,xlspan,y0); 
yOl = y 1 (xn 1 
yO 1(1,1) = Th; 

[x2,y2] = ode45('hothx2stl\x2span,y01);% Temperature at Th 

y02 = [y2(xn2,:) yin(8)]'; % Add the 8th element (i.e. H2a) for the ist stack 

[x3,y3] = ode45(’stackl',x3span,y02); 

y03 = y3(xn3,l:7); % only the first 7 elements are used in the hot heat exchanger 

[x4,y4] = ode45('ambhx2str,x4span,y03);% Target Tc at the 1st ambient heat exchanger 
y04 = y4(xn2,:); 

[x5,y5] = ode45('duct2stl’,x5span,y04); 
y05 = y5(xnl,:); 
y05(l,l) = Tc; 

[x6,y6] = ode45('ambhx2st2\x6span,y05);% Temperature at Tc 

y06 = [y6(xn2,:) yin(8) yin(9)]'; % Add the 8th element (i.e. H2) for the 2nd stack 

[x7,y7] = ode45('stack2’,x7span,y06); 

y07 = y7(xn3,l:7); % only the first 7 elements are used in the hot heat exchanger 

[x8,y8] = ode45('hothx2st2',x8span,y07); % Target Th at 2nd hot heat exchanger 
% Combine results of the integration 

% This is the result at the end of the 2nd hot heat exchanger 
yend = [y8(xn2,:) yin(8) yin(9)]'; 

% Compute the target values % 

fvec(l,l) = y03(l,l)-Tc; % Find the difference of the targeted cold temperature at ambhx2stl 
fvec(6,l) = yend(l,l)-Th;% Find the difference of the targeted hot temperature at hothx2st2 
% Find the differences of the target pre, pim, Ure and Uim % 
fvec(2:5,l)= yin(2:5,l)-yend(2:5,l); 

yend = yend(:); % Show results of the integration at the end 

% End of funcv2st.m % 



172 



APPENDIX D. PROPERTIES OF AF 45 GLASS 



Mechanical properties: 

Density at 20 °C (68 °F): p = 2.72 g/cm 3 

Young’s modulus: c = 66.0 kN/mm 2 
Torsional modulus: E = 26.7 kN/mm 2 

Poisson’s ratio: y= 0.235 
Electrical properties: 

Dielectric constant (1 MHz): E r = 6.2 
Dielectric loss factor (1 MHz): tan 8 = 9 x 10 4 

Temperature for a specific electrical resistivity of 10 8 Q cm: T Kl00 = 610 °C (1130 °F) 
Thermal properties: 

Viscosity and corresponding temperatures 



Viscosity 


Temperature 


Designation 


log T) (d Pas) 


in °C 


(°F) 




14.5 


627 


(1161) 


Strain point 


13.0 


663 


(1225) 


Annealing point 


7.60 


883 


(1621) 


Softening point 



Transformation temperature: T g = 662 °C (1224 °F) 

Coefficient of mean linear thermal expansion: (X 2 o- 3 oo = 4.5 x 10' 6 K' 1 
Optical properties: 

Refractive indices at 20°C (68°F): n e (A = 546 nm) = 1.5275, n d (A. = 546 nm) = 1.5255 
Abbe value: v 0 = 62.2% 

Luminous transmittance (glass thickness 1.1 mm): x VD65 = 91.2% 
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APPENDIX E.l CONSTRUCTION DRAWINGS FOR THE PRIME MOVER 

STACK HOLDER 



Horizontal bar 




0 . 048 +/- 0 . 002 ” 



0.248 +/- 
0 . 002 ” 




0 . 064 +/- 0 . 002 ” 

0 . 120 +/- 0 . 002 ” 

0 . 064 +/- 0 . 002 ” 



screw) 
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APPENDIX E.2 CONSTRUCTION DRAWINGS FOR THE AMBIENT 



HEAT EXCHANGER 



Vertical bar 
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APPENDIX E.2(CONTINUED) CONSTRUCTION DRAWINGS FOR THE 

AMBIENT HEAT EXCHANGER 

Horizontal bar 



0.148+/- 

0 . 002 ” 




0.047+/-0.002” 




0.074+/- 

0 . 002 ” 



drill-through 

(for 0-80 flat head screw) 



Central bar 



thickness : 0.200”+0.002” 



0.050+/-0.002” 







(- - 

1 .758+/-0.002” [ 

f ' 




drill-threaded (for 0-80 flat head screw) 
both ends 



10 slits, 0.025” deep, 0.028” wide 
with seperation of 0.160” (center to center) 
from edge to the first slit is 0.159” 
both sides 



0 . 100 ”+/- 0.002 
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APPENDIX E.3 CONSTRUCTION DRAWINGS FOR THE HOT HEAT 



EXCHANGER 




0.320+/-0.002” 



0 . 100 ”+/-- _ 
0 . 002 ” | 




(1) . 9 threaded holes of spacing 0.205”(center to center) for 0-80 3/16” fine screw 

(2) . 8 threaded holes of spacing 0.205” (center to center) for 0-80 3/16” fine screw 
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APPENDIX F. GRAPHS OF MATLAB AND MEASURED RESULTS 



This section which is divided into two subsections contains graphs of the results of 
the MATLAB program and measurements. First, results of the constricted annular prime 
mover are presented. Next, some predictied results for the two stack annular prime mover 
are provided. 
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APPENDIX F.l RESULTS OF THE CONSTRICTED ANNULAR PRIME 



MOVER 



Low mode, Resonance frequency vs. DeltaT , Rs = 0.7 




Figure F. 1-1 Comparisons of the measured and calculated resonance frequencies 
of the low mode for the constricted prime mover with the porosity of 0.7. 
Method A is the conformation transformation and Method B is the higher order 
mode. 
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High mode, Resonance frequency vs. DeltaT , Rs = 0.7 




Figure F.l-2 Comparisons of the measured and calculated resonance frequencies 
of the high mode for the constricted prime mover with the porosity of 0.7. 
Method A is the conformation transformation and Method B is the higher order 
mode. 
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Low mode, Resonance frequency vs. DeltaT , Rs = 0.3 




Figure F.l-3 Comparisons of the measured and calculated resonance frequencies 
of the low mode for the constricted prime mover with the porosity of 0.3. 
Method A is the conformation transformation and Method B is the higher order 
mode. 
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High mode, Resonance frequency vs. DeltaT , Rs = 0.3 




Figure F.l-4 Comparisons of the measured and calculated resonance frequencies 
of the high mode for the constricted prime mover with the porosity of 0.3. 
Method A is the conformation transformation and Method B is the higher order 
mode. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 0 Kelvin ; freq = 41 7.9+2. 149i ; Rs = 0.7 




Figure F.l-5 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 77.03 Kelvin ; freq = 420.5+1 .801 i ; Rs = 0.7 




Figure F.l-6 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=77 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 231 Kelvin ; freq = 429.1+0.5389i ; Rs = 0.7 




Figure F.l-7 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 0 Kelvin ; freq = 438.7+5.992i ; Rs = 0.7 




Figure F.l-8 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 77.03 Kelvin ; freq = 445+6.1 53i ; Rs = 0.7 




Figure F.l-9 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=77 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 231 Kelvin ; freq = 460+6.746i ; Rs = 0.7 




Figure F.l-10 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end corrections ; dT = 228.2 Kelvin ; freq = 307.4-1. 165i ; Rs = 0.1 




Figure F. 1-1 1 Mode shape of the low mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT =228 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. Symbol * represents the mode shape above onset, the 
measured frequency is 312 Hz. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 0.2 Kelvin ; freq = 361. 7+2.61 9i ; Rs = 0.3 




Figure F. 1-12 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 71 .73 Kelvin ; freq = 362.6+1 .847i ; Rs = 0.3 




Figure F.l-13 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=72 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end effect 2 ; dT = 226.5 Kelvin ; freq = 369.9+0.2682i ; Rs = 0.3 




Figure F.l-14 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 0.2 Kelvin ; freq = 466.9+6.538i ; Rs = 0.3 




Figure F.l-15 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 71.73 Kelvin ; freq = 468.7+6.1 94i ; Rs = 0.3 




Figure F.l-16 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=72 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



High mode with end effect 2 ; dT = 226.5 Kelvin ; freq = 479.6+5.239i ; Rs = 0.3 




Figure F.l-17 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 

calculated results are based on the higher order mode method. The frequency is 
the calculated value. 
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Normalized pressure 



Low mode with end effect 3 ; dT = 0 Kelvin ; freq = 295.6+3.873i ; Rs = 0.1 




Figure F. 1-18 Mode shape of the low mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



Low mode with end effect 3 ; dT = 228.2 Kelvin ; freq = 306.8-0.6823i ; Rs = 0.1 




Figure F. 1-19 Mode shape of the low mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=228 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



High mode with end effect 3 ; dT = 0 Kelvin ; freq = 468.8+9.638i ; Rs = 0.1 




Figure F.l-20 Mode shape of the high mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



High mode with end effect 3 ; dT = 228.2 Kelvin ; freq = 488.9+8.246i ; Rs = 0.1 




Figure F. 1-21 Mode shape of the high mode of the constricted prime mover 
(Rs=0.1) when the driver is located 45° from the stack and AT=228 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



Low mode with end effect 3 ; dT = 0.2 Kelvin ; freq = 361. 3+2.881 i ; Rs = 0.3 




Figure F.l-22 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



Low mode with end effect 3 ; dT = 226.5 Kelvin ; freq = 369.6+0.4252i ; Rs = 0.3 




Figure F.l-23 Mode shape of the low mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



High mode with end effect 3 ; dT = 0.2 Kelvin ; freq = 461.9+9.625i ; Rs = 0.3 




Figure F.l-24 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



High mode with end effect 3 ; dT = 226.5 Kelvin ; freq = 475.5+8.55i ; Rs = 0.3 




Figure F.l-25 Mode shape of the high mode of the constricted prime mover 
(Rs=0.3) when the driver is located 45° from the stack and AT=227 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



Low mode with end effect 3 ; dT = 0 Kelvin ; freq = 417.8+2.217i ; Rs = 0.7 




Figure F.l-26 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



Low mode with end effect 3 ; dT = 231 Kelvin ; freq = 428.8+0.4826i ; Rs = 0.7 




Figure F.l-27 Mode shape of the low mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



High mode with end effect 3 ; dT = 0 Kelvin ; freq = 434.2+8.779i ; Rs = 0.7 




Figure F.l-28 Mode shape of the high mode of the constricted prime mover 
(Rs-0.7) when the driver is located 45° from the stack and AT=0 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



High mode with end effect 3 ; dT = 231 Kelvin ; freq = 456.3+9. 995i ; Rs = 0.7 




Figure F.l-29 Mode shape of the high mode of the constricted prime mover 
(Rs=0.7) when the driver is located 45° from the stack and AT=231 K. The 

calculated results are based on the higher order mode method and the end 
corrections of the stack are included. The frequency is the calculated value. 
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Normalized pressure 



APPENDIX F.2 RESULTS OF THE TWO STACK ANNULAR 

PRIME MOVER 



Mode shape for of the two stack prime mover, dT = 25K; freq = 425+7.421 i 




Figure F.2-1 The calculated mode shape of the high mode of the two stack 
annular prime mover at AT = 25 K. Case 1 : The duct between the two hot heat 
exchangers is held at T a . 
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Normalized pressure 



Mode shape for of the two stack prime mover, dT = -5.466e-05K; freq = 424+7.234i 




Figure F.2-2 The calculated mode shape of the low mode of the two stack 
annular prime mover at AT = 0 K. Case 1: The duct between the two hot heat 
exchangers is held at T a . 
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